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The  contribution  of  unresolved  sources  to  the  diffuse  gamma-ray  background  could  induce  anisotropies 
in  this  emission  on  small  angular  scales.  We  analyze  the  angular  power  spectrum  of  the  diffuse  emission 
measured  by  the  Fermi  Large  Area  Telescope  at  Galactic  latitudes  \b\  >  30°  in  four  energy  bins  spanning 
1-50  GeV.  At  multipoles  (  >  155,  corresponding  to  angular  scales  S  2°,  angular  power  above  the  photon 
noise  level  is  detected  at  >99.99%  confidence  level  in  the  1-2  GeV,  2-5  GeV,  and  5-10  GeV  energy  bins, 
and  at  >99%  confidence  level  at  10-50  GeV.  Within  each  energy  bin  the  measured  angular  power  takes 
approximately  the  same  value  at  all  multipoles  €  >  155,  suggesting  that  it  originates  from  the  contribution 
of  one  or  more  unclustered  source  populations.  The  amplitude  of  the  angular  power  normalized  to  the 
mean  intensity  in  each  energy  bin  is  consistent  with  a  constant  value  at  all  energies,  CP/(/)2  =  9.05  ± 

0.84  X  10-6  sr,  while  the  energy  dependence  of  CP  is  consistent  with  the  anisotropy  arising  from  one  or 
more  source  populations  with  power-law  photon  spectra  with  spectral  index  Ts  =  2.40  ±  0.07.  We  discuss 
the  implications  of  the  measured  angular  power  for  gamma-ray  source  populations  that  may  provide  a 
contribution  to  the  diffuse  gamma-ray  background. 
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I.  INTRODUCTION 

The  origin  of  the  all-sky  diffuse  gamma-ray  emission 
remains  one  of  the  outstanding  questions  in  high-energy 
astrophysics.  First  detected  by  OSO-3  [1],  the  isotropic 
gamma-ray  background  (IGRB)  was  subsequently  mea¬ 
sured  by  SAS-2  [2],  the  Energetic  Gamma  Ray 
Experiment  Telescope  [3,4],  and  most  recently  by  the 
Large  Area  Telescope  (LAT)  on  board  the  Fermi 
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Gamma-Ray  Space  Telescope  (Fermi)  [5].  The  term  IGRB 
is  used  to  refer  to  the  observed  diffuse  gamma-ray  emission 
which  appears  isotropic  on  large  angular  scales  but  may 
contain  anisotropies  on  small  angular  scales.  The  IGRB 
describes  the  collective  emission  of  unresolved  members 
of  extragalactic  source  classes  and  Galactic  source  classes 
that  contribute  to  the  observed  emission  at  high  latitudes, 
and  gamma-ray  photons  resulting  from  the  interactions  of 
ultra-high-energy  cosmic  rays  with  intergalactic  photon 
fields  [6]. 

Confirmed  gamma-ray  source  populations  with 
resolved  members  are  guaranteed  to  contribute  to  the 
IGRB  at  some  level  via  the  emission  from  fainter,  un¬ 
resolved  members  of  those  source  classes.  In  the 
Energetic  Gamma  Ray  Experiment  Telescope  era  the 
possibility  that  blazars  are  the  dominant  contributor  to 
the  IGRB  intensity  was  extensively  studied  (e.g.,  [7-9]), 
however  the  level  of  the  blazar  contribution  remains 
uncertain,  with  recent  results  suggesting  different 
energy-dependent  contributions  from  blazars,  which 
amount  to  as  little  as  —15%  or  as  much  as  —100%  of 
the  Fermi-measured  IGRB  intensity,  depending  on  the 
energy  [10-13].  Star-forming  galaxies  [14]  and  gamma- 
ray  millisecond  pulsars  [15]  may  also  provide  a  signifi¬ 
cant  contribution  to  the  IGRB  at  some  energies. 
However,  substantial  uncertainties  in  the  properties  of 
even  confirmed  source  populations  present  a  challenge 
to  estimating  the  amount  of  emission  attributable  to 
each  source  class,  and  currently  the  possibility  that  the 
IGRB  includes  an  appreciable  contribution  from  un¬ 
known  or  unconfirmed  gamma-ray  sources,  such  as 
dark  matter  annihilation  or  decay  (e.g.,  [16-22]),  cannot 
be  excluded. 

The  Fermi-measured  IGRB  energy  spectrum  is  rela¬ 
tively  featureless,  following  a  simple  power  law  to  good 
approximation  over  a  large  energy  range  (—200  MeV  to 
—  100  GeV)  [5],  As  a  result,  identifying  the  contribu¬ 
tions  from  individual  components  based  on  spectral 
information  alone  is  difficult.  However,  in  addition  to 
the  energy  spectrum  and  average  intensity,  the  IGRB 
contains  angular  information  in  the  form  of  fluctuations 
on  small  angular  scales  [23].  The  statistical  properties 
of  these  small-scale  anisotropies  may  be  used  to  infer 
the  presence  of  emission  from  unresolved  source 
populations. 

If  some  component  of  the  IGRB  emission  originates 
from  an  unresolved  source  population,  rather  than  from  a 
perfectly  isotropic,  smooth  source  distribution,  the  dif¬ 
fuse  emission  will  contain  fluctuations  on  small  angular 
scales  due  to  the  varying  number  density  of  sources 
in  different  sky  directions.  Unlike  the  Poisson  fluctua¬ 
tions  between  pixels  in  a  map  of  a  truly  isotropic 
source  distribution  (which  we  shall  call  “photon  noise”), 
which  are  due  to  finite  event  statistics,  the  fluctuations 
from  an  unresolved  source  population  are  inherent 
in  the  source  distribution  and  will  not  decrease  in 
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amplitude  even  in  the  limit  of  infinite  statistics.  Hence, 
with  sufficient  statistics,  these  fluctuations  could  be 
detected  above  those  expected  from  the  photon  noise, 
and  could  be  used  to  understand  the  origin  of  the  diffuse 
emission. 

The  angular  power  spectrum  of  the  emission  provides  a 
metric  for  characterizing  the  intensity  fluctuations.  For  a 
source  population  modeled  with  a  specific  spatial  and 
luminosity  distribution,  the  angular  power  spectrum  can 
be  predicted  and  compared  to  the  measured  angular  power 
spectrum;  in  this  way  an  anisotropy  measurement  has  the 
potential  to  constrain  the  properties  of  source  populations. 
Other  approaches  to  using  anisotropy  information  in  the 
IGRB  have  also  been  considered.  For  example,  the  1 -point 
probability  distribution  function  (PDF),  i.e.  the  distribution 
of  the  number  of  counts  per  pixel,  is  an  alternative  metric  to 
characterize  the  fluctuations  [13,24,25].  In  addition,  cross- 
correlating  the  gamma-ray  sky  with  galaxy  catalogs  or  the 
cosmic  microwave  background  can  be  used  to  constrain  the 
origin  of  the  emission  [26]. 

In  recent  years  theoretical  studies  have  predicted  the 
angular  power  spectrum  of  the  gamma-ray  emission 
from  several  known  and  proposed  source  classes. 
Established  astrophysical  source  populations  such  as 
blazars  [27-29],  star-forming  galaxies  [30],  and 
Galactic  millisecond  pulsars  [31]  have  been  considered 
as  possible  contributors  to  the  anisotropy  of  the  IGRB. 
In  addition,  it  has  been  shown  that  the  annihilation  or 
decay  of  dark  matter  in  Galactic  subhalos  [32-34]  and 
extragalactic  structures  [23,27,29,34-39],  may  generate 
an  anisotropy  signal  in  diffuse  gamma-ray  emission. 
Interestingly,  the  predicted  angular  power  spectra  of 
these  gamma-ray  source  classes  in  the  multipole  range 
of  L  —  100-500  are  in  most  cases  fairly  constant  in 
multipole  (except  for  dark  matter  annihilation  and  decay 
signals,  e.g.,  [23,27,38]),  although  the  amplitude  of  the 
predicted  anisotropy  varies  between  source  classes.  This 
multipole-independent  signal  arises  from  the  Poisson 
term  in  the  angular  power  spectrum,  which  describes 
the  anisotropy  from  an  unclustered  collection  of  point 
sources.  The  multipole  independence  of  the  predicted 
angular  power  spectra  therefore  indicates  that  the  ex¬ 
pected  degree  of  intrinsic  clustering  of  these  gamma-ray 
source  populations  has  a  subdominant  effect  on  the 
angular  power  spectra  in  this  multipole  range.  The  an¬ 
gular  power  spectra  of  dark  matter  annihilation  and 
decay  signals  are  predicted  to  be  smooth  and  relatively 
featureless,  with  the  angular  power  generally  falling  off 
more  quickly  with  multipole  than  Poisson  angular 
power. 

In  this  work  we  present  a  measurement  of  the  angular 
power  spectrum  of  the  high-latitude  emission  detected  by 
the  Fermi  LAT,  using  —22  months  of  data.  The  data  were 
processed  with  the  FERMI  SCIENCE  TOOLS  [40],  and  binned 
into  maps  covering  several  energy  ranges.  Regions  of  the 
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sky  heavily  contaminated  by  Galactic  diffuse  emission 
and  known  point  sources  were  masked,  and  then  angular 
power  spectra  were  calculated  on  the  masked  sky  for  each 
energy  bin  using  the  HEALPIX  package  [41],  described 
in  [42], 

To  understand  the  impact  of  the  instrument  response 
on  the  measured  angular  power  spectrum,  several  tai¬ 
lored  validation  studies  were  performed  for  this  analysis. 
The  robustness  of  the  anisotropy  analysis  pipeline  was 
tested  using  a  source  model  with  known  anisotropy 
properties  that  was  simulated  to  include  the  effects  of  the 
instrument  response  and  processed  with  the  same  analy¬ 
sis  pipeline  as  the  data.  The  data  processing  was  cross¬ 
checked  to  exclude  the  presence  of  anisotropies  created 
by  systematics  in  the  instrument  exposure  calculation  by 
using  an  event-shuffling  technique  (as  used  in  [43])  that 
does  not  rely  on  the  Monte-Carlo-based  exposure  calcu¬ 
lation  implemented  in  the  SCIENCE  TOOLS.  In  addition, 
validation  studies  were  performed  to  characterize  the 
impact  of  foreground  contamination,  masking,  and  inac¬ 
curacies  in  the  assumed  point  spread  function  (PSF). 

We  use  a  set  of  simulated  models  of  the  gamma-ray  sky 
as  a  reference,  and  compare  the  angular  power  spectrum 
measured  for  the  data  to  that  of  the  models  to  identify  any 
significant  differences  in  anisotropy  properties.  Finally,  we 
compare  the  predicted  anisotropy  for  several  confirmed 
and  proposed  gamma-ray  source  populations  to  the  mea¬ 
sured  angular  power  spectrum  of  the  data. 

The  data  selection  and  map-making  procedure  are 
described  in  Sec.  II,  and  the  angular  power  spectrum 
calculation  is  outlined  in  Sec.  III.  The  event-shuffling 
technique  is  presented  in  Sec.  IV,  and  the  details  of  the 
models  simulated  to  compare  with  the  data  are  given  in 
Sec.  V.  The  results  of  the  angular  power  spectrum  mea¬ 
surement  and  the  validation  studies  are  presented  in 
Sec.  VI.  The  energy  dependence  of  the  anisotropy  is  dis¬ 
cussed  in  Sec.  VII,  and  the  implications  of  the  results  for 
specific  source  populations  are  examined  in  Sec.  VIII.  The 
conclusions  are  summarized  in  Sec.  IX. 

II.  DATA  SELECTION  AND  PROCESSING 

The  Fermi  LAT  is  designed  to  operate  primarily  as 
a  survey  instrument,  featuring  both  a  wide  field  of  view 
(—2.4  sr)  and  a  large  effective  area  (£7000  cm2  for 
normally-incident  photons  above  1  GeV).  The  telescope 
is  equipped  with  a  4  X  4  array  of  modules,  each  consisting 
of  a  precision  tracker  and  calorimeter,  covered  by  an  anti- 
coincidence  detector  that  allows  for  rejection  of  charged 
particle  events.  Full  details  of  the  instrument,  including 
technical  information  about  the  detector,  onboard  and 
ground  data  processing,  and  mission-oriented  support,  are 
given  in  [44], 

We  selected  data  taken  from  the  beginning  of  scientific 
operations  in  early- August  2008  through  early- June  2010, 
encompassing  over  56.6  Ms  of  live  time  [45].  We  selected 


only  “diffuse”  class  [44]  events  to  ensure  that  the  events 
are  photons  with  high  probability,  and  restricted  our  analy¬ 
sis  to  the  energy  range  1-50  GeV  where  the  PSF  of  the 
LAT  is  small  enough  to  allow  for  sufficient  sensitivity  to 
anisotropies  at  small  angular  scales.  The  upper  limit  of 
50  GeV  was  chosen  because  the  small  photon  statistics 
above  this  energy  severely  limit  the  sensitivity  of  the 
analysis  at  the  high  multipoles  of  interest.  The  data  and 
simulations  were  analyzed  with  the  LAT  analysis  software 
SCIENCE  TOOLS  version  v9rl5p4  using  the  standard  P6_V3 
LAT  instrument  response  functions  (IRFs).  Detailed  docu¬ 
mentation  of  the  SCIENCE  TOOLS  is  given  in  [46]. 

In  order  to  both  promote  near  uniform  sky  exposure 
and  to  limit  contamination  from  gamma  rays  originating 
in  Earth’s  atmosphere,  the  tool  gtrnktime  was  used  to 
remove  data  taken  during  any  time  period  when  the  LAT 
rocked  to  an  angle  exceeding  52°  with  respect  to  the 
zenith,  and  during  any  time  period  when  the  LAT  was 
not  in  survey  mode.  Beginning  in  its  second  year  of 
operation  (September  2009),  Fermi  has  been  operating 
in  survey  mode  with  a  large  rocking  angle  of  50°,  in 
contrast  to  the  35°  rocking  angle  used  during  the  first 
year  of  operation.  The  rocking-angle  cut  is  used  to  limit 
the  amount  of  contamination  from  gamma  rays  produced 
in  cosmic -ray  interactions  in  the  upper  atmosphere  by 
using  only  data  taken  when  the  Earth’s  limb  was  outside 
of  the  field  of  view  (the  Earth’s  limb  has  zenith  angle 
—  113°).  However,  due  to  the  LAT’s  large  field  of  view, 
some  Earth-limb  gamma  rays  may  be  observed  even 
when  the  rocking  angle  constraint  is  not  exceeded,  thus 
the  gtselect  tool  was  also  used  to  remove  each  individual 
event  with  a  zenith  angle  exceeding  105°.  We  note  that 
all  events  in  the  data  set  were  detected  while  the  Fermi 
spacecraft  was  outside  of  the  South  Atlantic  Anomaly 
region  in  which  the  cosmic-ray  fluxes  at  the  altitude  of 
Fermi  are  significantly  enhanced. 

In  order  to  balance  the  need  for  a  large  effective  area 
with  the  need  for  high  angular  resolution,  the  LAT  uses  a 
combination  of  thin  tracker  regions  near  the  front  of  the 
instrument  and  thicker  tracker  regions  in  the  back  of  the 
detector.  While  the  effective  area  of  each  region  is  compa¬ 
rable,  the  width  of  the  PSF  for  events  detected  in  the  front 
trackers  is  approximately  half  that  of  events  detected  in  the 
back  of  the  instrument.  For  a  measurement  of  the  angular 
power  at  high  multipoles,  it  is  thus  necessary  to  differ¬ 
entiate  between  photons  observed  in  the  front  and  back 
trackers  of  the  Fermi  LAT.  In  this  study,  we  processed 
front-  and  back-converting  events  separately,  using  the 
gtselect  tool  to  isolate  each  set  of  events  and  calculating 
the  exposure  maps  independently.  The  P6_V3_DIFFUSE: 
FRONT  and  P6JV3„DIFFUSE:BACK  IRFs  were  used  to  analyze 
the  corresponding  sets  of  events. 

Taking  the  selection  cuts  into  account,  the  integrated  live 
time  was  calculated  using  gtltcube.  We  chose  a  pixel  size 
of  0.125°,  which  produces  a  HEALPIX  map  with  resolution 
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parameter  Aside  =  512.  At  this  resolution,  the  suppression 
of  angular  power  from  the  pixel  window  function  is  sub¬ 
dominant  with  respect  to  the  suppression  from  the  LAT 
PSF.  We  adopted  an  angular  step  size  cos(0)  =  0.025  in 
order  to  finely  grid  the  exposure  map  for  different  gamma- 
ray  arrival  directions  in  instrument  coordinates.  The 
exposure  was  then  calculated  using  gtexpcube  with 
the  same  pixel  size,  for  42  logarithmic  energy  bins  span¬ 
ning  1.04-50.0  GeV.  These  finely-gridded  energy  bins 
were  then  summed  to  build  maps  covering  four  larger 
energy  bins,  as  described  in  Sec.  Ill  A.  Using  the 
GaRDiAn  package  [47],  both  the  photon  counts  and  expo¬ 
sure  maps  were  converted  into  HE ALPIX- format  maps  with 

^sid  e  =  512. 

III.  ANGULAR  POWER  SPECTRUM 
CALCULATION 

We  consider  the  angular  power  spectrum  Q  of  an  in¬ 
tensity  map  1(f)  where  ip  denotes  the  sky  direction.  The 
angular  power  spectrum  is  given  by  the  coefficients 
Q  =  (\afm\2),  with  the  a%m  determined  by  expanding  the 
map  in  spherical  harmonics, 

/(<A)  =  Y^a,mYM).  (1) 

im 

The  intensity  angular  power  spectrum  indicates  the  dimen¬ 
sionful  size  of  intensity  fluctuations  and  can  be  compared 
with  predictions  for  source  classes  whose  collective  inten¬ 
sity  is  known  or  assumed  (as  in,  e.g.,  [31]).  The  intensity 
angular  power  spectrum  of  a  single  source  class  is  not  in 
general  independent  of  energy  due  to  the  energy  depen¬ 
dence  of  the  mean  map  intensity  (I). 

We  can  also  construct  the  dimensionless  fluctuation 
angular  power  spectrum  by  dividing  the  intensity  angular 
power  spectrum  Q  of  a  map  by  the  mean  sky  intensity 
(outside  of  the  mask,  for  a  masked  sky  map)  squared  (I)2. 
The  fluctuation  angular  power  spectrum  characterizes  the 
angular  distribution  of  the  emission  independent  of  the 
intensity  normalization.  Its  amplitude  for  a  single  source 
class  is  the  same  in  all  energy  bins  if  all  members  of  the 
source  class  share  the  same  observed  energy  spectrum, 
since  this  results  in  the  angular  distribution  of  the  collec¬ 
tive  emission  being  independent  of  energy.  Energy  depen¬ 
dence  in  the  fluctuation  angular  power  due  to  variation  of 
the  energy  spectra  between  individual  members  of  the 
population  is  discussed  in  Sec.  VII. 

A.  Energy  dependence 

We  calculate  the  angular  power  spectrum  of  the  data  and 
simulated  models  in  four  energy  bins.  Using  multiple 
energy  bins  increases  sensitivity  to  source  populations 
that  contribute  significantly  to  the  anisotropy  in  a  limited 
energy  range,  and  may  also  aid  in  the  interpretation  of  a 
measurement  in  terms  of  a  detection  of  or  constraints  on 
specific  source  populations  [39,48].  In  addition,  the  detec¬ 


tion  of  an  energy  dependence  in  the  fluctuation  angular 
power  spectrum  of  the  total  emission  (the  anisotropy  en¬ 
ergy  spectrum)  may  be  used  to  infer  the  presence  of 
multiple  contributing  source  classes  [49].  In  the  case  that 
a  single  source  population  dominates  the  anisotropy  over  a 
given  energy  range,  the  energy  dependence  of  the  intensity 
angular  power  spectrum  can  indicate  the  energy  spectrum 
of  that  contributor. 

Since  the  LAT’s  angular  resolution  and  the  photon  sta¬ 
tistics  depend  strongly  on  energy,  the  sensitivity  of  the 
analysis  is  also  energy-dependent:  at  low  energies  the 
LAT’s  PSF  broadens,  resulting  in  reduced  sensitivity  to 
small-scale  anisotropies,  while  at  high  energies  the  mea¬ 
surement  uncertainties  are  dominated  by  low  statistics. 
We  calculate  angular  power  spectra  in  the  energy  bins 
1.04-1.99  GeV,  1.99-5.00  GeV,  5.00-10.4  GeV,  and 
10.4-50.0  GeV.  The  map  for  each  energy  bin  for  the 
angular  power  spectrum  analysis  was  created  by  summing 
the  corresponding  maps  produced  in  finely-gridded  energy 
bins,  as  described  in  Sec.  II. 

B.  Angular  power  spectrum  of  a  masked  sky 

The  focus  of  this  work  is  to  search  for  anisotropies  on 
small  angular  scales  from  unresolved  source  populations, 
hence  the  regions  of  the  sky  used  in  this  analysis  were 
selected  to  minimize  the  contribution  of  the  Galactic  dif¬ 
fuse  emission  from  cosmic-ray  interactions  and  the  emis¬ 
sion  from  known  sources.  A  mask  excluding  Galactic 
latitudes  \b\  <  30°  and  a  2°  angular  radius  around  each 
source  in  the  11-month  Fermi  LAT  catalog  (1FGL)  [50] 
was  applied  prior  to  performing  the  angular  power  spec¬ 
trum  calculations  in  all  energy  bins.  The  fraction  of  the  sky 
outside  of  this  mask  is  /sky  =  0.325.  The  2°  angular  radius 
for  the  source  masking  approximately  corresponds  to  the 
95%  containment  angle  for  events  at  normal  incidence  at 
1  GeV  (front/back  average  for  P6_V3  IRFs);  the  contain¬ 
ment  angle  decreases  with  increasing  energy.  The  effect  of 
the  mask  on  the  angular  power  spectra  is  discussed  below 
and  in  Sec.  VI F,  and  the  impact  of  variations  in  the  latitude 
cut  is  assessed  in  Sec.  VI E.  An  all-sky  intensity  map  of  the 
data  in  each  energy  bin  is  shown  in  Fig.  1,  both  with  and 
without  applying  the  default  mask. 

The  angular  power  spectra  of  the  masked  maps  were 
calculated  using  HEALPIX,  after  first  removing  the  mono¬ 
pole  and  dipole  terms.  To  approximately  correct  for  the 
power  suppression  due  to  masking,  the  raw  angular  power 
spectra  output  by  HEALPIX  were  divided  by  the  fraction  of 
the  sky  outside  the  mask,  /sky.  This  correction  is  valid  at 
multipoles  greater  than  —10,  where  the  power  spectrum  of 
the  signal  varies  much  more  slowly  than  the  window 
function,  as  detailed  below. 

When  a  fraction  of  the  sky  is  masked,  the  measured 
spherical  harmonics  coefficients  are  related  to  the  true, 
underlying  spherical  harmonics  coefficients,  off,  via  a 
matrix  multiplication,  a?m  =  where 
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DATA  (P6_V3  diffuse),  1. 0-2.0  GeV 


-7.0  -4.0  Log  (Intensity  [cm  2  s  1  sr  ‘J) 


DATA  (P6_V3  diffuse),  2.0 -5.0  GeV 


DATA  (P6_V3  diffuse),  2. 0-5.0  GeV 


-4.0  Log  (Intensity  [cm  2  s  1  sr  ']) 


-4.0  Log  (Intensity  [cm  2  s  1  sr  ']) 


i  -4.0  Log  (Intensity  [cm  2  s  1  sr  ']) 


DATA  (P6_V3  diffuse),  5.0-10.4  GeV 


DATA  (P6_V3  diffuse),  5.0-10.4  GeV 


—4.0  Log  (Intensity  [ci 


— 4.0  Log  (Intensity  [ci 


DATA  (P6_V3  diffuse),  10.4-50.0  GeV 


DATA  (P6_V3  diffuse),  10.4-50.0  GeV 


FIG.  1  (color  online).  All-sky  intensity  maps  of  the  data  in  the  four  energy  bins  used  in  this  analysis,  in  Galactic  coordinates;  the  map 
projection  is  Mollweide.  The  data  shown  are  the  average  of  the  maps  of  the  front-  and  back-converting  events,  and  are  shown 
unmasked  (left  panels)  and  with  the  default  mask  applied  (right  panels).  The  mask  excludes  Galactic  latitudes  \b\  <  30°  and  a  2° 
angular  radius  around  each  source  in  the  1FGL  catalog.  The  map  images  shown  have  been  downgraded  in  resolution  to  iVside  =  128  to 
improve  the  visual  quality  of  the  images;  however,  the  analysis  was  performed  on  the  higher  resolution  maps  as  described  in  the  text. 
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Wu'mmf  is  so-called  coupling  matrix  given  by 
=  /nobs  d2nY*m(n)Y(,m,{ n),  where  the  integral  is 
done  only  on  the  unmasked  sky  whose  solid  angle  is  (lobs. 
Then,  HEALPIX  returns  a  raw  angular  power  spectrum, 
C™w  =  (2£  +  ir‘X,„|a€m|2,  whose  ensemble  average  is 
related  to  the  true  power  spectrum,  C/Llc ,  as 

<qaw>  =  2€'+  (2) 

mm' 

Now,  for  a  given  mask,  one  may  calculate 
and  estimate  C‘ILIe  by  inverting  this  equation.  This  approach 
is  called  the  MASTER  algorithm  [51],  and  has  been  shown 
to  yield  unbiased  estimates  of  C‘ILle .  However,  while  this  is 
an  unbiased  estimator,  it  is  not  necessarily  a  minimum- 
variance  one.  In  particular,  when  the  coupling  matrix  is 
nearly  singular  because  of,  e.g.,  an  excessive  amount  of 
mask  or  a  complex  morphology  of  mask,  this  estimator 
amplifies  noise.  We  observed  this  amplification  of  noise 
when  applying  the  MASTER  algorithm  to  our  data  set. 
Therefore,  we  decided  to  use  an  approximate,  but  less 
noisy  alternative.  It  is  easy  to  show  (see,  e.g.,  Eq.  (A3)  of 
[52])  that,  when  Xmm'lIkfFmm'l2  peaks  sharply  at  l  =  E 
and  C[rLlc  varies  much  more  slowly  than  the  width  of  this 
peak,  the  above  equation  can  be  approximated  as 

<cr>  «  qr  ^  =  que/sky.  (3) 

This  approximation  eliminates  the  need  for  a  matrix  inver¬ 
sion.  We  have  verified  that  this  method  yields  an  unbiased 
result  with  substantially  smaller  noise  than  the  MASTER 
algorithm  at  €  >  10.  We  adopt  this  method  throughout  this 
paper. 


by  Ref.  [53]  (see  Eq.  [20]  of  that  work)  for  a  three- 
dimensional  density  field,  but  the  same  is  true  for  a  two- 
dimensional  field,  as  we  are  dealing  with  here.  We  have 
verified  this  using  numerical  simulations. 

In  this  paper,  although  we  use  maps  at  Aside  =  512  (for 
which  the  maximum  multipole  is  €max  =  1024),  we  restrict 
the  analysis  to  Q  up  to  f max  —  500  where  we  have  a 
reasonable  signal-to-noise  ratio.  For  these  multipoles  the 
effect  of  the  pixel  window  function  is  negligible,  and  thus 
we  shall  simplify  our  analysis  pipeline  by  not  applying  the 
pixel  window  correction  to  the  observed  power  spectrum 
[54].  Therefore,  our  signal  power  spectrum  estimator  is 
given  by 


^signal 


£raw  / 


//sky  Q\ 


(Wt5™)2 


(4) 


where  CN  =  <A/,pix>(lM2ix>/npix  is  the  photon  noise 
term,  with  Nyp ix,  Apix,  and  Opjx  the  number  of  observed 
events,  the  exposure,  and  the  solid  angle,  respectively,  of 
each  pixel,  and  the  averaging  is  done  over  the  unmasked 
pixels.  We  approximate  the  photon  noise  term  by  CN  = 
( I)24Trfsky/Ny ,  with  Ny  denoting  the  total  number  of 
observed  events  outside  the  mask.  This  approximation  is 
accurate  at  the  percent  level.  Note  that  while  qaw  is  always 
non-negative,  it  is  possible  for  our  estimator  for  the  signal 
power  spectrum  C/gnal  to  be  negative  due  to  the  subtraction 
of  the  noise  term. 

The  beam  window  function  in  multipole  space  associ¬ 
ated  with  the  full  non-Gaussian  PSF  is  given  by 


W/eam(E)  =  2tt 


dcos0P({cos{d))PSF{d\E), 


(5) 


C.  Window  functions 

The  angular  power  spectrum  calculated  from  a  map  is 
affected  by  the  PSF  of  the  instrument  and  the  pixelization 
of  the  map,  encoded  in  the  beam  window  function  W/J,eam 
and  the  pixel  window  function  W|1X  respectively,  both  of 
which  can  lead  to  a  multipole-dependent  suppression  of 
angular  power  that  becomes  stronger  at  larger  multipoles. 
Depending  upon  whether  the  power  spectrum  originates 
from  signal  or  noise,  corrections  for  the  beam  and  pixel 
window  functions  must  be  applied  to  the  measurement 
differently.  For  our  application,  we  must  not  apply  any 
corrections  to  the  photon  shot  noise  (Poisson  noise)  term, 
while  we  must  apply  both  the  beam  and  pixel  window 
function  corrections  to  the  signal  term  from,  e.g.,  unre¬ 
solved  sources.  While  it  is  obvious  why  one  must  not  apply 
the  beam  window  correction  to  the  photon  noise  term,  it 
may  not  be  so  obvious  why  one  must  also  not  apply  the 
pixel  window  correction  to  that  same  term.  In  fact,  this 
statement  is  correct  only  for  the  shot  noise,  if  the  data  are 
pixelized  by  the  nearest-grid  assignment  (which  we  have 
adopted  for  our  pixelizing  scheme).  This  has  been  shown 


where  Pf{cos{9))  are  the  Legendre  polynomials  and 
PSF(0;  E)  is  the  energy-dependent  PSF  for  a  given  set  of 
IRFs,  with  9  denoting  the  angular  distance  in  the  distribu¬ 
tion  function.  The  PSF  used  corresponds  to  the  average  for 
the  actual  pointing  and  live  time  history  of  the  LAT  and 
over  the  off-axis  angle,  as  given  by  the  gtpsf  tool.  We 
calculate  the  beam  window  functions  for  both  the  front- 
and  back-converting  events. 

The  PSF  of  the  LAT,  and  consequently  the  beam  window 
function,  varies  substantially  over  the  energy  range  used  in 
this  analysis,  and  also  non-negligibly  within  each  energy 
bin.  We  treated  this  energy  dependence  by  calculating  an 
average  window  function  <VPr^eam(J5’,))  for  each  energy  bin 
Ej,  weighted  by  the  intensity  spectrum  of  the  events  in  each 
bin, 

(W]?eam (£■,•))  =  — [E™‘  dEWfDE)^,  (6) 

4in  JE^t  dt 

where  7bin  =  /f™"'  dE(dN / dE)  and  Emm/l  and  Emax  i  are 
the  lower  and  upper  edges  of  each  energy  bin.  The  differ¬ 
ential  intensity  dN/dE  outside  the  mask  in  each  map  for 
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the  finely-gridded  energy  bins  described  in  Sec.  II  was  used 
to  approximate  the  energy  spectrum  for  this  calculation. 


D.  Measurement  uncertainties 

The  ler  statistical  uncertainty  crC(  on  the  measured 
angular  power  spectrum  coefficients  C“gnal  is  given  by  [55] 


(2€  +  l)/skyA€ 


(cfnal 


C, 


(W>?eam)2 


> 


(7) 


where  A€  is  the  width  of  the  multipole  bin  (for  binned 
data). 

After  implementing  the  corrections  for  masking  and  for 
the  beam  window  function  to  estimate  the  signal  angular 
power  spectrum  via  Eq.  (4),  the  coefficients  Cj'8’1'11  were 
binned  in  multipole  with  A£  =  50  and  averaged  in  each 
multipole  bin,  weighted  by  the  measurement  uncertainties, 


<Q> 


IQ/< 

TT747 


(B) 


with  Cf  =  Cfnal  as  calculated  by  Eq.  (4)  and  ctq  given  by 
Eq.  (7)  with  A€  =  1  and  VT>?eam  =  (lT^eam(£())  for  the 
corresponding  energy  bin  E,.  As  expected,  we  find  that 
the  statistical  measurement  uncertainties  calculated  at  the 
linear  center  of  each  multipole  bin  via  Eq.  (7)  with  A€  = 
50  agree  well  with  the  scatter  within  each  multipole  bin. 
The  value  of  Q  at  multipoles  2  <  €  <  4  was  found  in  most 
cases  to  be  anomalously  large  [56],  indicating  the  presence 
of  strong  correlations  on  very  large  angular  scales,  such  as 
those  that  could  be  induced  by  the  shape  of  the  mask  and 
by  contamination  from  Galactic  diffuse  emission.  To  avoid 
biasing  the  value  of  the  average  C“gnal  in  the  first  bin  by  the 
values  at  these  low  multipoles,  the  multipole  bins  begin  at 
€  =  5. 

Finally,  the  angular  power  spectra  of  the  front-  and  back- 
converting  events  were  combined  by  weighted  averaging, 
weighting  by  the  measurement  uncertainty  on  each  data 
point.  Because  of  the  larger  PSF  associated  with  back- 
converting  events,  the  measurement  errors  on  the  angular 
power  spectra  of  the  back-converting  data  set  tend  to  be 
larger  than  those  of  the  front-converting  data  set,  particu¬ 
larly  at  low  energies  and  high  multipoles  where  the  sup¬ 
pression  of  the  raw  angular  power  due  to  the  beam  window 
function  is  much  stronger  for  the  back-converting  data  set. 
The  difference  between  the  measurement  uncertainties 
associated  with  the  front  and  back  data  sets  is  less  promi¬ 
nent  at  higher  energies. 


IV.  EVENT-SHUFFLING  TECHNIQUE 

One  way  to  search  for  anisotropies  is  to  first  calculate  the 
flux  of  particles  from  each  direction  in  the  sky  (equal  to  the 
number  of  detected  events  from  some  direction  divided  by 


the  exposure  in  the  same  direction),  and  then  examine  its 
directional  distribution.  The  flux  calculation,  which  re¬ 
quires  knowledge  of  the  exposure,  depends  on  the  effective 
area  of  the  detector  and  the  accumulated  observation  live 
time. 

The  effective  area,  calculated  from  a  Monte  Carlo  simu¬ 
lation  of  the  instrument,  could  suffer  from  systematic 
errors,  such  as  miscalculations  of  the  dependence  of  the 
effective  area  on  the  instrument  coordinates  (off-axis  angle 
and  azimuthal  angle).  Naturally,  any  systematic  errors 
involved  in  the  calculation  of  the  exposure  will  propagate 
to  the  flux,  possibly  affecting  its  directional  distribution.  If 
the  magnitude  of  these  systematic  errors  is  comparable  to 
or  larger  than  the  statistical  power  of  the  available  data  set, 
their  effects  on  the  angular  distribution  of  the  flux  might 
masquerade  as  a  real  detectable  anisotropy.  For  this  reason, 
we  cross-check  our  results  using  an  alternative  method  to 
construct  an  exposure  map  that  does  not  rely  on  the  Monte- 
Carlo-based  calculation  of  the  exposure  implemented  in 
the  SCIENCE  TOOLS. 

The  starting  point  of  this  method  is  the  construction  of  a 
sky  map  that  shows  how  an  isotropic  sky  would  look  as 
seen  by  the  Fermi  LAT.  This  sky  map,  hereafter  called  the 
“no-anisotropy  sky  map,”  is  directly  proportional  to  the 
exposure  map. 

One  method  of  generating  a  no-anisotropy  map  is  to 
randomize  the  reconstructed  directions  of  the  detected 
events  (as  in  [43]).  In  the  case  that  the  angular  distribution 
of  the  flux  is  perfectly  isotropic,  a  time-independent  inten¬ 
sity  should  be  detected  when  looking  in  any  given  detector 
direction.  Possible  time  variation  of  the  intensity  would  be 
due  only  to  changes  in  the  operating  conditions  of  the 
instrument.  A  set  of  isotropic  events  can  be  built  by  ran¬ 
domly  coupling  the  times  and  the  directions  of  real  events 
in  local  instrument  coordinates.  The  randomization  in  this 
analysis  was  performed  by  exchanging  the  direction  of  a 
given  real  event  in  the  LAT  frame  with  the  direction  of 
another  event  selected  randomly  from  the  data  set  with 
uniform  probability.  Using  this  information,  the  sky  direc¬ 
tion  is  reevaluated  for  the  two  events.  By  construction,  the 
randomized  data  set  preserves  the  exposure,  the  energy  and 
angular  (with  respect  to  the  LAT  reference  frame)  distri¬ 
butions,  and  also  accounts  for  the  detector  dead  times. 

As  already  discussed  in  Sec.  II,  for  this  analysis  a  cut  of 
52°  on  the  rocking  angle  was  applied  to  limit  possible 
photon  contamination  from  the  Earth’s  albedo.  For  the 
shuffling  technique,  the  analysis  was  performed  with  a 
reduced  field  of  view  of  the  instrument,  namely,  the  events 
used  were  selected  to  have  an  off-axis  angle  less  than  50°. 
In  this  way,  events  with  zenith  angle  exceeding  102°  were 
removed.  This  selection  cut  avoids  introducing  asymme¬ 
tries  in  the  exposure  across  the  held  of  view  due  to  cutting 
events  based  on  zenith  angle. 

The  randomization  was  performed  using  the  masked  sky 
map  described  in  Sec.  Ill,  so  that  only  real  events  with  sky 
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coordinates  outside  the  masks  were  used,  and  the  reeval¬ 
uated  sky  direction  for  each  event  was  required  to  be  in  the 
unmasked  region  of  the  sky.  This  randomization  process 
was  repeated  20000  times,  separately  for  the  front-  and 
back-converting  events,  each  time  producing  a  shuffled  sky 
map  that  is  compatible  with  an  isotropic  source  distribu¬ 
tion.  The  final  no-anisotropy  sky  map  for  each  energy  bin 
was  produced  by  taking  the  average  of  these  20000 
shuffled  sky  maps.  For  the  available  event  statistics,  aver¬ 
aging  20000  shuffled  maps  was  reasonably  effective  at 
reducing  the  Poisson  noise  associated  with  the  average 
number  of  events  per  pixel.  To  reduce  the  number  of 
required  shuffled  maps  by  increasing  the  average  number 
of  events  per  pixel,  the  shuffled  maps  were  constructed  at 
slightly  lower  resolution  (N^  =  256)  than  was  used  in  the 
default  analysis.  When  analyzing  the  anisotropy  with  these 
exposure  maps  from  the  shuffling  technique,  count  maps  at 
Aside  =  256  were  used  to  construct  the  intensity  maps.  A 
no-anisotropy  sky  map  is  shown  in  Fig.  2.  This  sky  map 
does  not  appear  entirely  uniform  because  the  sky  was  not 
observed  with  uniform  exposure. 

Although  the  no-anisotropy  sky  map  is  directly  propor¬ 
tional  to  the  exposure  map,  this  method  does  not  allow  us 
to  determine  the  absolute  level  of  the  exposure.  We  there¬ 
fore  constructed  intensity  maps  (with  arbitrary  normaliza¬ 
tion)  by  dividing  the  real  data  counts  maps  in  each  energy 
bin  by  the  no-anisotropy  map  for  that  energy  bin,  after  first 
smoothing  the  no-anisotropy  map  with  a  Gaussian  beam 
with  er  =  1°  to  reduce  the  pixel-to-pixel  fluctuations  due 
to  the  finite  number  of  events  available  to  use  in  the 
randomization.  This  smoothing  beam  size  removes  noise 
in  the  no-anisotropy  sky  map  above  €  ~  200,  and  was 
chosen  because  we  focus  our  search  for  anisotropies  in 
that  multipole  range.  Angular  power  spectra  were  then 
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FIG.  2  (color  online).  No-anisotropy  sky  map  created  by  sum¬ 
ming  20000  shuffled  maps  using  front-  and  back-converting 
events  with  E  >  1  GeV,  binned  into  a  HEALPIX  map  with  /Vside  = 
256.  The  map  projection  is  Hammer-Aitoff.  The  features  in  the 
no-anisotropy  sky  map  result  from  the  fact  that  the  sky  was  not 
observed  with  uniform  exposure. 
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calculated  from  these  intensity  maps  as  in  Sec.  III. 
Because  of  the  arbitrary  normalization  of  these  intensity 
maps,  we  calculate  only  fluctuation  angular  power  spectra 
of  the  data  when  using  the  exposure  map  produced  by  this 
shuffling  technique. 

V.  SIMULATED  MODELS 

Detailed  Monte  Carlo  simulations  of  Fermi  LAT  all-sky 
observations  were  performed  to  provide  a  reference  against 
which  to  compare  the  results  obtained  for  the  real  data  set. 
The  simulations  were  produced  using  the  gtobssim  tool, 
which  simulates  observations  with  the  LAT  of  an  input 
source  model.  The  gtobssim  tool  generates  simulated  pho¬ 
ton  events  for  an  assumed  spacecraft  pointing  and  live-time 
history,  and  a  given  set  of  IRFs.  The  P6_V3_DIFFUSE  IRFs 
and  the  actual  spacecraft  pointing  and  live-time  history 
matching  the  observational  time  interval  of  the  data  were 
used  to  generate  the  simulated  data  sets. 

Two  models  of  the  gamma-ray  sky  were  simulated.  Each 
model  is  the  sum  of  three  components: 

(1)  GAL — A  model  of  the  Galactic  diffuse  emission. 

(2)  CAT — The  sources  in  the  1 1 -month  catalog  (1FGL) 
[50]. 

(3)  ISO — An  isotropic  background. 

Both  models  include  the  same  CAT  and  ISO  compo¬ 
nents,  and  differ  only  in  the  choice  of  the  model  for  the 
GAL  component.  GAL  describes  both  the  spatial  distribu¬ 
tion  and  the  energy  spectrum  of  the  Galactic  diffuse  emis¬ 
sion.  The  GAL  component  for  the  reference  sky  model 
used  in  this  analysis  (hereafter,  MODEL)  is  the  recom¬ 
mended  Galactic  diffuse  model  for  LAT  data  analysis, 
GLL_IEM_V02.FIT  [57],  which  has  an  angular  resolution  of 
0.5°.  This  model  was  used  to  obtain  the  1LGL  catalog;  a 
detailed  description  can  be  found  in  Ref.  [58]. 

An  alternate  sky  model  (ALT  MODEL)  was  simulated 
for  comparison,  in  order  to  test  the  possible  impact  of 
variations  in  the  Galactic  diffuse  model.  This  model  is 
internal  to  the  LAT  collaboration,  and  was  built  using  the 
same  method  as  GLL_IEM_V02.FIT,  but  differs  primarily  in 
the  following  ways:  (i)  this  model  was  constructed  using 
21  months  of  Fermi  LAT  observations,  while 
GLL_IEM_V02.FIT  was  based  on  9  months  of  data,  and 
(ii)  additional  large-scale  structures,  such  as  the  Lermi 
bubbles  [59],  are  included  in  the  model  through  the  use 
of  simple  templates. 

The  sources  in  CAT  were  simulated  with  energy  spectra 
approximated  by  single  power  laws,  and  with  the  locations, 
average  integral  fluxes,  and  photon  spectral  indices  as 
reported  in  the  1LGL  catalog.  All  1451  sources  were 
included  in  the  simulation.  ISO  represents  the  sum 
of  the  Lermi-measured  IGRB  and  an  additional  isotropic 
component  presumably  due  to  unrejected  charged 
particles;  for  this  component  the  spectrum  template 
ISOTROPIC_IEM_V02.TXT  was  used. 
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For  both  the  MODEL  and  the  ALT  MODEL,  the  sum  of 
the  three  simulated  components  results  in  a  description  of 
the  gamma-ray  sky  that  closely  approximates  the  angular- 
dependent  intensity  and  energy  spectrum  of  the  all-sky 
emission  measured  by  the  Fermi  LAT.  Although  the  simu¬ 
lated  models  may  not  accurately  reproduce  some  large- 
scale  structures,  e.g.,  Loop  I  [60]  and  the  Fermi  bubbles, 
these  features  are  not  expected  to  induce  anisotropies  on 
the  small  angular  scales  on  which  we  focus  in  this  work. 

VI.  RESULTS 

In  this  section  we  present  the  measured  angular  power 
spectra  of  the  data,  followed  by  the  results  of  validation 
studies  which  examine  the  effect  of  variations  in  the  de¬ 
fault  analysis  parameters,  and  by  a  comparison  of  the 
results  for  the  data  with  those  for  simulated  models.  We 
summarize  the  main  results  of  the  angular  power  spectrum 
measurements  of  the  data  and  of  key  validation  studies 
described  in  this  section  in  Table  I. 

Unless  otherwise  noted,  the  results  are  shown  for  data 
and  models  with  the  angular  power  spectra  calculated  after 
applying  the  default  source  mask  which  excludes  sources 
in  the  1FGL  catalog  and  Galactic  latitudes  |Z?|<30°. 
Because  of  the  arbitrary  normalization  of  the  intensity 
maps  calculated  using  the  exposure  map  from  shuffling, 
we  show  fluctuation  angular  power  spectra  for  this  data  set. 
Intensity  angular  power  spectra  are  presented  for  all  other 
data  sets. 

In  the  figures  we  show  our  signal  angular  power  spec¬ 
trum  estimator  Csf's":l  [Eq.  (4)],  which  represents  the  signal 
after  correcting  for  the  power  suppression  due  to  masking, 


subtracting  the  photon  noise,  and  correcting  for  the  beam 
window  function.  A  measurement  that  is  inconsistent  with 
zero  thus  indicates  the  presence  of  signal  angular  power. 
The  Cfnal  shown  is  the  weighted  average  of  this  quantity 
for  the  maps  of  front  and  back  events.  The  fluctuation 
angular  power  spectra  C“gnal/(/)2  were  calculated  by  di¬ 
viding  Cfnal  of  the  front  and  back  events  by  their  respec¬ 
tive  (I)2,  and  then  averaging  the  angular  power  spectra.  For 
conciseness,  in  the  figure  labels  Q  =  C[aw//sky  is  the  raw 
angular  power  spectrum  output  by  HEALPIX  corrected  for 
the  effects  of  masking.  The  error  bars  on  points  indicate  the 
1  er  statistical  uncertainty  in  the  measurement  in  each  mul¬ 
tipole  bin  as  calculated  by  Eq.  (7)  with  Af  =  50  and  with 
the  bins  beginning  at  f  =  5.  The  binned  data  points  are 
located  at  the  linear  center  of  each  multipole  bin. 

A.  Angular  power  spectrum  of  the  data 

We  now  present  the  results  of  the  angular  power  spec¬ 
trum  analysis  of  the  data.  We  measure  the  angular  power 
spectrum  of  the  data  after  applying  the  default  latitude  cut 
and  source  mask,  and  refer  to  this  as  our  default  data 
analysis  (DATA).  We  also  measure  the  angular  power 
spectrum  of  the  data  using  the  same  masking  and  analysis 
pipeline  after  performing  Galactic-foreground  cleaning, 
described  below,  and  refer  to  this  as  the  cleaned  data 
analysis  (DATA: CLEANED).  These  two  measurements 
constitute  our  main  results  for  the  data,  and  so  we  discuss 
the  energy  dependence  of  the  measured  angular  power 
(Sec.  VII)  and  present  constraints  on  specific  source  pop¬ 
ulations  (Sec.  VIII)  for  the  results  of  both  the  default  and 
cleaned  data  analyses. 


TABLE  I.  Best- fit  values  of  the  angular  power  CP  and  fluctuation  angular  power  CP/(/)2  in  each  energy  bin  over  the  multipole  range 
155  <  €  <  504.  Results  are  shown  for  the  data  processed  with  the  default  analysis  pipeline,  the  foreground-cleaned  data,  the  data 
analyzed  with  the  2FGL  source  mask,  and  the  default  simulated  model.  Significance  indicates  the  measured  angular  power  expressed 
in  units  of  the  measurement  uncertainty  cr;  the  measurement  uncertainties  can  be  taken  to  be  Gaussian. 


Emm  [GeV] 

£max  [GeV] 

CP  [(cm  2  s  1  sr  *)2  srl 

Significance 

CP/(/)2  [10“6  sr] 

Significance 

DATA 

1.04 

1.99 

7.39  ±  1.14  X  10“18 

6.5  a 

10.2  ±  1.6 

6.5  a 

1.99 

5.00 

1.57  ±  0.22  X  10“ 18 

1.2a 

8.35  ±  1.17 

1.1a 

5.00 

10.4 

1.06  ±  0.26  X  10“19 

4.1  a 

9.83  ±  2.42 

4. 1  a 

10.4 

50.0 

2.44  ±  0.92  X  10“20 

2.1a 

8.00  ±  3.37 

2.4  a 

DATA:CLEANED 

1.04 

1.99 

4.62  ±  1.11  X  10“18 

4.2  a 

6.38  ±  1.53 

4.2  a 

1.99 

5.00 

1.30  ±  0.22  X  10“ 18 

6.0  a 

6.90  ±  1.16 

5.9  a 

5.00 

10.4 

8.45  ±  2.46  X  10“20 

3.4  a 

8.37  ±  2.41 

3.5  a 

10.4 

50.0 

2.11  ±  0.86  X  10“20 

2.4  a 

7.27  ±  3.36 

2.2  a 

DATA:2FGL 

1.04 

1.99 

5.18  ±  1.17  X  10“18 

4.4  a 

7.23  ±  1.61 

4.5  a 

1.99 

5.00 

1.21  ±  0.28  X  10“18 

5.3cr 

6.49  ±  1.22 

5.3  a 

5.00 

10.4 

8.38  ±  2.72  X  10“20 

3.1  a 

7.67  ±  2.54 

3.0  a 

10.4 

50.0 

8.00  ±  9.57  X  10“21 

0.8  a 

2.28  ±  3.52 

0.6  a 

MODEL 

1.04 

1.99 

1.89  ±  1.08  X  10“18 

0.1a 

2.53  ±  1.47 

1.1a 

1.99 

5.00 

1.92  ±  2.10  X  10“19 

0.9  a 

0.99  ±  1.12 

0.9  a 

5.00 

10.4 

3.41  ±  2.60  X  10“20 

1.3  a 

3.04  ±  2.34 

1.3  a 

10.4 

50.0 

0.62  ±  9.63  X  10“21 

0.1  a 

0.24  ±  3.02 

0.  lcr 
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To  minimize  the  impact  of  Galactic  foregrounds  we  have 
employed  a  large  latitude  cut.  However,  Galactic  diffuse 
emission  extends  to  very  high  latitudes  and  may  not  exhibit 
a  strong  gradient  with  latitude,  and  it  is  thus  important  to 
investigate  to  what  extent  our  data  set  may  be  contaminated 
by  a  residual  Galactic  contribution.  For  this  purpose  we 
attempt  to  reduce  the  Galactic  diffuse  contribution  to  the 
high-latitude  emission  by  subtracting  a  model  of  the 
Galactic  foregrounds  from  the  data,  and  then  calculating 
the  angular  power  spectra  of  the  residual  maps.  For  the 
angular  power  spectrum  analysis  of  the  residual  maps 
(cleaned  data)  we  note  that  the  noise  term  CN  is  calculated 
from  the  original  (uncleaned)  map,  since  subtracting  the 
model  from  the  data  does  not  reduce  the  photon  noise  level. 

In  the  following  we  use  the  recommended  Galactic  dif¬ 
fuse  model  GLL_IEM_V 02. FIT,  which  is  also  the  default  GAL 
model  that  we  simulate,  as  described  in  Sec.  V.  To  tailor  the 
model  to  the  high-latitude  sky  regions  considered  in  this 
work,  the  normalization  of  the  model  was  adjusted  by  refit¬ 
ting  the  model  to  the  data  only  in  the  regions  outside  the 
latitude  mask.  For  the  fit  we  used  GaRDiAn  which  con¬ 
volves  the  model  with  the  instrument  response  (effective 
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area  and  PSF).  The  normalization  obtained  in  this  way  is, 
however,  very  close  to  the  nominal  one,  within  a  few  percent. 

We  present  the  angular  power  spectra  of  the  data  before 
and  after  Galactic-foreground  cleaning  in  Fig.  3;  expanded 
versions  of  the  angular  power  spectra  for  the  1-2  GeV  and 
2-5  GeV  bins  focusing  on  the  high-multipole  data  are 
shown  in  Fig.  4.  In  both  analyses,  angular  power  at  €  > 
155  is  measured  in  the  data  in  all  energy  bins  considered, 
and  the  angular  power  spectra  for  the  default  and  cleaned 
data  are  in  good  agreement  in  this  multipole  range.  In  the 
default  data,  the  large  increase  in  angular  power  at  i  <  155 
in  the  two  energy  bins  spanning  1-5  GeV  is  likely  due  to 
contamination  from  the  Galactic  diffuse  emission  which 
features  correlations  on  large  angular  scales,  but  may  also 
be  attributable  in  part  to  the  effects  of  the  source  mask  (see 
Sec.  VI F). 

At  €  >  155  the  measured  angular  power  does  not  exhibit 
a  clear  scale  dependence  in  any  energy  bin.  The  results  of 
fitting  the  unbinned  signal  angular  power  spectrum  estima¬ 
tor  for  155  <  i  <  504  in  each  energy  bin  to  a  power  law 
(-.signal  rf-  (€/€q)”  with  fn  =  155  are  given  in  Table  II  for  the 
default  data  analysis.  In  each  energy  bin,  the  angular  power 
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FIG.  3  (color  online).  Comparison  of  intensity  angular  power  spectra  of  the  data  and  Galactic-foreground-cleaned  data.  For  €  >  155 
the  measured  power  at  all  energies  is  approximately  constant  in  multipole,  suggesting  that  it  originates  from  one  or  more  unclustered 
source  populations.  The  large  increase  in  angular  power  in  the  default  data  at  f  <  155  in  the  1-2  and  2-5  GeV  bins  is  likely  attributable 
largely  to  contamination  from  Galactic  diffuse  emission.  In  these  two  energy  bins,  foreground  cleaning  primarily  reduces  angular 
power  at  €  <  155,  with  the  most  significant  reductions  at  (  <  105.  At  energies  greater  than  5  GeV  the  effect  of  foreground  cleaning  is 
small  for  (  >  55.  Expanded  versions  of  the  top  panels  are  shown  in  Fig.  4. 
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FIG.  4  (color  online).  Expanded  versions  of  top  panels  of  Fig.  3,  focusing  on  the  high-multipole  angular  power. 


spectrum  for  155  <  €  <  504  is  consistent  with  a  Poisson 
spectrum  (constant  in  multipole,  i.e.,  n  =  0  falls  within  the 
95%  confidence  level  (CL)  range  of  the  best-fit  power-law 
index),  as  expected  for  the  angular  power  spectrum  of  one 
or  more  unclustered  source  populations.  However,  we 
emphasize  that  the  uncertainty  in  the  scale  dependence  is 
appreciable,  particularly  for  the  10-50  GeV  bin. 

In  light  of  the  scale  independence  of  the  angular  power 
at  f  >  155,  we  associate  the  signal  in  this  multipole  range 
with  a  Poisson  angular  power  spectrum  and  determine  the 
best-fit  constant  value  of  the  angular  power  CP  and  the 
fluctuation  angular  power  CP/(/)2  over  155  <  €  <  504  in 
each  energy  bin,  by  weighted  averaging  of  the  unbinned 
measurements.  These  results  for  the  default  and  cleaned 
data  are  summarized  in  Table  I,  along  with  the  results 
obtained  for  the  data  using  an  updated  source  catalog  to 
define  the  source  mask  and  for  a  simulated  model,  which 
will  be  discussed  in  Sec.  VI G  and  VI H  respectively. 

We  note  that  the  associated  measurement  uncertainties 
can  be  taken  to  be  Gaussian,  in  which  case  the  reported 
significance  quantifies  the  probability  of  the  measured 
angular  power  to  have  resulted  by  chance  from  a  truly 
uniform  background.  We  consider  a  3cr  or  greater  detection 
of  angular  power  (CP)  in  a  single  energy  bin  to  be  statis¬ 
tically  significant.  For  the  default  data,  the  best-fit  values  of 
CP  indicate  significant  detections  of  angular  power  in  the 
1-2,  2-5,  and  5-10  GeV  bins  (6.5cr,  1.2a,  and  4. lcr,  re¬ 
spectively),  while  in  the  10-50  GeV  bin  the  best-fit  CP 
represents  a.  2.1  a  measurement  of  angular  power.  We  further 
note  that  the  best-fit  value  of  the  fluctuation  angular  power 
over  all  four  energy  bins  (see  Sec.  VII)  yields  a  detection 
with  greater  than  10cr  significance  for  the  default  data. 

For  the  1-2  GeV  and  2-5  GeV  energy  bands  the  cleaning 
procedure  results  in  a  significant  decrease  in  the  angular 
power  at  low  multipoles  (€  <  105),  and  a  smaller  reduction 
at  higher  multipoles.  However,  the  decrease  is  small  for 
€  >  155,  and  angular  power  is  still  measured  at  all  ener¬ 
gies,  at  slightly  lower  significances  (see  Table  I).  We 
emphasize  that  the  detections  in  the  three  energy  bins 
spanning  1-10  GeV  remain  statistically  significant,  and 


the  best-fit  fluctuation  angular  power  over  all  energy  bins 
is  detected  at  greater  than  8cr  significance.  For  energies 
above  5  GeV  the  foreground  cleaning  does  not  strongly 
affect  the  measured  angular  power  spectrum  for  €  >  55.  At 
all  energies  the  decrease  in  angular  power  at  low  multi¬ 
poles  can  be  attributed  to  the  reduction  of  Galactic  fore¬ 
grounds  which  feature  strong  correlations  on  large  angular 
scales.  We  conclude  that  contamination  of  the  data  by 
Galactic  diffuse  emission  does  not  have  a  substantial  im¬ 
pact  on  our  results  at  the  multipoles  of  interest  (€  >  155). 
This  conclusion  is  in  agreement  with  that  of  Ref.  [39], 
which  found  that  the  Galactic  foregrounds  have  a  rapidly 
declining  angular  power  spectrum  above  €  ~  100. 

To  further  study  the  expected  angular  power  spectrum  of 
Galactic  foregrounds,  we  analyzed  the  angular  power  spec¬ 
trum  of  the  E(B-V)  emission  map  of  Ref.  [61]  (hereafter 
SFD  map),  which  is  proportional  to  the  column  density  of 
the  interstellar  dust,  after  masking  \b\  <  30°  as  in  our 
default  analysis.  The  SFD  map  is  a  good  tracer  of  the 
Galactic  interstellar  medium  (ISM)  away  from  the 
Galactic  plane,  the  spatial  structure  of  which  should  be 
reflected  in  the  diffuse  gamma-ray  emission  produced  by 
interactions  of  cosmic  rays  with  the  ISM.  It  has  an  angular 
resolution  of  6  arcminutes,  much  smaller  than  the  intrinsic 
resolution  of  the  GAL  model  map  (  ~  0.5°),  and  smaller 
than  the  map  resolution  used  in  this  study,  and  so  it 
accurately  represents  the  small-scale  structure  of  the  ISM 
on  the  angular  scales  accessible  to  this  analysis.  We  found 

TABLE  II.  Multipole  dependence  of  intensity  angular  power 
in  the  data  (default  analysis)  for  155  <  (.  <  504  in  each  energy 
bin.  The  best-fit  power-law  index  n  in  each  energy  bin  is  given 
with  the  associated  \ 2  Per  degree  of  freedom  (d.o.f.)  of  the  fit. 


F 

^min 

F 

^max 

n 

*2/d.o.f. 

1.04 

1.99 

-1.33  ±  0.78 

0.38 

1.99 

5.00 

-0.07  ±  0.45 

0.43 

5.00 

10.4 

-0.79  ±  0.76 

0.37 

10.4 

50.0 

-1.54  ±  1.15 

0.39 

083007-12 


ANISOTROPIES  IN  THE  DIFFUSE  GAMMA-RAY  . . . 


PHYSICAL  REVIEW  D  85,  083007  (2012) 


that  the  SFD  map  produces  an  angular  power  spectrum 
with  a  slightly  harder  slope  than  the  default  GAL  model, 
and  consequently  features  more  angular  power  at  high 
multipoles.  However,  like  the  GAL  model,  the  SFD  map 
angular  power  spectrum  falls  off  quickly  with  multipole 
compared  to  a  Poisson  spectrum,  and  the  amplitude  of  the 
SFD  map  angular  power  is  below  that  measured  in  the  data 
for  €  S:  100.  This  further  reinforces  the  conclusion  that 
Galactic  foreground  contamination  cannot  explain  the  ob¬ 
served  high-multipole  angular  power  in  the  data. 


B.  Validation  with  a  simulated  point  source  population 

To  ensure  that  our  analysis  procedure  accurately  recov¬ 
ers  an  input  angular  power  spectrum,  and,  in  particular,  that 
the  result  is  not  biased  by  instrumental  effects,  we  compare 
the  angular  power  spectrum  calculated  for  a  simulated 
point  source  population  with  the  theoretical  prediction 
for  that  population.  It  is  straightforward  to  calculate  the 
expected  angular  power  spectrum  of  unclustered  point 
sources,  once  a  flux  distribution  function,  dN/dS  (in  units 
of  cm2  s  sr  1 ),  and  a  source  detection  flux  threshold,  Sc 
(in  units  of  cm-2  s_l),  are  provided.  The  angular  power 
spectrum  of  an  unclustered  point  source  population  is  the 
Poisson  component  of  the  angular  power  CP,  which  takes 
the  same  value  at  all  multipoles  and  is  given  by 

fs,  dN 

CP=  /  dSS2  — .  (9) 

Jo  dS 


For  our  source  population  model  we  adopt  the  best-fit 
flux  distribution  for  the  high-latitude  Fermi  sources,  re¬ 
ported  in  [10],  which  describes  dN/dS  with  a  broken 
power-law  model: 


dN 

~dS 


AS~P',  S  >  Sh, 
AS~pi+plS~P\  S<Sb, 


(10) 


which  contains  four  free  parameters,  A,  /S x ,  02,  and  Sb.  For 
this  form  of  dN/dS,  the  source  power  spectrum  can  be 
found  analytically  (for  Sc  >  Sb): 


CP  =  A 


S3~* 

3-Pi 


0i  ~  02  (SbV-Pil 
3  -  02  \SC)  j 


(11) 


A  fit  for  the  simulated  source  population  for  1 .04- 
10.4  GeV  yields  A  =  (1.90  ±  0.48)  X  10“13(180/7r)2, 
0i  =  2.213  ±  0.073,  02  =  1.533  ±  0.007,  and  Sb  = 
1.41  X  10  9  cm  2  s  1 .  Note  that  the  errors  are  correlated. 
Figure  5  shows  the  predicted  CP  for  this  source  model  for 
threshold  fluxes  in  the  range  Sc  =  0.8  X  1 0  7  —  4.2  X 
10~7  cm"2s-1.  The  error  bars  are  calculated  from  the 
full  covariance  matrix  of  the  above  parameters.  Although 
we  have  used  zero  as  the  lower  limit  of  the  integral  in 
Eq.  (9),  using  the  actual  lower  limit  of  the  flux  distribution 
adopted  for  the  simulated  population  results  in  a  negligible 
difference  in  the  predicted  CP. 


FIG.  5.  Predicted  amplitude  of  the  source  angular  power  spec¬ 
trum,  CP  [see  Eq.  (9)],  for  energies  of  1.04-10.4  GeV  as  a 
function  of  a  source  detection  threshold  flux,  Sc. 


We  simulated  this  source  population  model  with  gtobs- 
sim  using  the  same  procedure  as  described  in  Sec.  V.  The 
simulated  population  comprises  nearly  20000  point 
sources  distributed  randomly  across  the  entire  sky,  with 
each  source’s  flux  drawn  from  the  flux  distribution  speci¬ 
fied  above.  The  photon  spectrum  of  each  source  is  modeled 
as  a  power  law  with  a  spectral  index  T  ( dN/dE  E~r) 
drawn  from  a  Gaussian  distribution  with  mean  of  2.40  and 
a  standard  deviation  of  0.28.  The  simulated  events  were 
processed  and  the  angular  power  spectrum  of  this  source 
model  calculated  using  the  same  procedure  as  was  used  for 
the  data  and  other  simulations  in  this  study,  except  that  the 
energy  range  of  the  map  was  chosen  to  be  1.04-10.4  GeV, 
and  no  mask  was  applied. 

The  fluxes  of  the  —20  000  simulated  sources  were  drawn 
from  a  flux  distribution  in  which  the  maximum  possible 
flux  ( E  >100  MeV)  that  could  be  assigned  to  a  source  was 
10-5  cm-2  s_1,  however  the  maximum  flux  of  any  source 
in  the  simulation,  which  represents  a  single  realization  of 
this  source  population,  was  —3  X  10-6  citT2  s-1.  We  take 
these  values  as  the  upper  and  lower  bound  on  the  source 
detection  threshold  flux  (E  >  100  MeV)  corresponding  to 
the  simulated  model,  since  we  do  not  impose  a  source 
detection  threshold  by  masking  or  otherwise  excluding 
simulated  sources  above  a  specific  threshold  flux.  A 
spectral  index  T  =  2.4  is  assumed  to  determine  the  thresh¬ 
old  fluxes  in  the  1.04-10.4  GeV  energy  band.  From  these 
threshold  fluxes  we  calculate  the  corresponding  upper  and 
lower  bound  on  the  predicted  CP  in  the  1.04-10.4  GeV 
energy  band. 

The  angular  power  spectrum  for  the  simulated  source 
population  calculated  via  the  analysis  pipeline  used  in  this 
study  is  presented  in  Fig.  6,  with  the  shaded  region  indicat¬ 
ing  the  predicted  range  of  CP  (the  mean  values  of  CP  at  the 
upper  and  lower  flux  threshold);  for  a  given  model  CP  is 
independent  of  multipole,  thus  we  expect  the  recovered 
angular  power  spectrum  to  be  independent  of  multipole 
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FIG.  6  (color  online).  Intensity  angular  power  spectrum  of 
a  simulated  observation  of  the  source  population  model, 
compared  with  the  theoretical  prediction  {shaded  band).  The 
angular  power  spectrum  of  the  simulated  population  is  in 
excellent  agreement  with  the  prediction  over  a  large  multipole 
range. 


with  amplitude  within  the  shaded  region.  The  angular 
power  spectrum  recovered  from  the  simulated  data  is  in 
excellent  agreement  with  the  prediction  up  to  multipoles  of 
€  ~  800.  Above  €  ~  800,  the  upturn  in  the  measured  an¬ 
gular  power  spectrum  is  likely  due  to  inaccuracies  in  the 
modeling  of  the  beam  window  function,  which  can  intro¬ 
duce  features  on  very  small  angular  scales.  In  the  remain¬ 
der  of  this  study,  we  present  results  only  for  the  multipole 
range  €  =  5  to  €  =  504. 

C.  Sensitivity  to  the  exposure  map  calculation 

To  investigate  the  possibility  that  potential  inaccuracies 
in  the  exposure  map  calculation  for  the  default  analysis 
might  generate  spurious  anisotropy  in  the  intensity  maps, 
we  compare  the  fluctuation  angular  power  spectra  of  the 
data  using  our  default  analysis  pipeline  with  the  results 
obtained  after  replacing  the  default  exposure  map  with  that 
generated  by  the  event-shuffling  technique  described  in 
Sec.  IV.  This  is  shown  in  Figs.  7  and  8.  In  these  two  figures 
only,  the  results  from  the  default  data  analysis  were  ob¬ 
tained  from  maps  of  HEALPIX  resolution  Vside  =  256  to 
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FIG.  7  (color  online).  Fluctuation  angular  power  spectra  Q/(/}2  calculated  using  the  default  analysis  pipeline  compared  with  those 
obtained  using  the  exposure  map  from  the  event  shuffling  technique  described  in  Sec.  IV.  Angular  power  is  measured  in  all  four  energy 
bins  by  both  analysis  methods.  The  lack  of  significant  differences  at  the  multipoles  of  interest  between  the  angular  power  spectra 
yielded  by  the  two  methods  demonstrates  that  any  inaccuracies  in  the  exposure  map  have  a  negligible  impact  on  the  measured  angular 
power  spectra.  Expanded  versions  of  the  top  panels  are  shown  in  Fig.  8. 
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FIG.  8  (color  online).  Expanded  versions  of  top  panels  of  Fig.  7. 


match  the  resolution  of  the  maps  using  the  exposure  de¬ 
termined  from  the  shuffling  technique.  All  other  results 
presented  in  this  study  were  obtained  from  /VSKle  =  512 
maps.  Because  of  the  reduced  map  resolution,  the  pixel 
window  function  has  a  small  effect  on  the  angular  power 
spectra  shown  in  Figs.  7  and  8,  however  it  affects  the 
results  of  the  default  analysis  and  the  analysis  using  the 
shuffled  exposure  map  in  the  same  way,  and  so  these  results 
can  still  be  directly  compared  for  the  purpose  of  checking 
the  effect  of  the  exposure  map  calculation. 

The  results  of  the  two  analysis  methods  are  in  good 
agreement  at  all  energies  and  multipoles  considered,  ex¬ 
cept  for  slight  deviations  at  €  <  55  for  1-5  GeV.  We 
caution  that  at  these  low  multipoles  the  measured  angular 
power  spectra  may  be  strongly  affected  by  the  mask,  which 
has  features  on  large  angular  scales.  The  slight  differences 
in  the  data  selection  cuts  for  the  analysis  using  the  expo¬ 
sure  map  from  the  shuffling  technique  compared  to  those 
for  the  default  data  analysis  could  lead  to  the  observed 


differences  in  the  low-multipole  angular  power  spectra. 
The  differences  could  also  result  from  systematics  in  the 
Monte-Carlo-based  exposure  calculation  implemented  in 
the  SCIENCE  TOOLS,  leading  to  inaccuracies  in  the  exposure 
map  which  vary  on  large  angular  scales.  As  we  do  not 
focus  on  the  low-multipole  angular  power  in  this  study,  we 
defer  a  full  investigation  of  this  issue  to  future  work.  The 
agreement  at  L  >  55  demonstrates  that  any  potential 
spatially-dependent  inaccuracies  in  the  SCIENCE  TOOLS 
exposure  calculation  have  a  negligible  impact  on  the  an¬ 
gular  power  spectra  in  the  multipole  range  of  interest.  In 
particular,  from  the  consistency  of  the  two  methods  we 
conclude  that  using  the  Monte-Carlo-based  exposure  cal¬ 
culation  does  not  induce  spurious  signal  anisotropy  in  our 
results. 

D.  Dependence  on  the  PSF  model 

We  examine  the  impact  of  variations  in  the  assumed  PSF 
on  the  results  of  the  analysis  by  comparing  the  beam 


FIG.  9  (color  online).  Comparison  of  the  beam  window  functions  for  the  P6_V3  and  P6_V1 1  IRFs;  the  P6_V3  IRFs  are  the  default  used 
in  this  analysis.  The  quantity  Wf ,  which  is  the  factor  by  which  the  angular  power  is  suppressed  due  to  the  finite  angular  resolution  of 
the  instrument,  is  shown  for  the  front-converting  ( left  panel)  and  back-converting  ( right  panel)  events,  evaluated  at  the  log-center  of 
each  energy  bin  used  in  this  analysis.  The  differences  between  the  of  these  two  IRFs  are  small  (  S  few  percent)  at  all  energies 
considered,  indicating  that  our  results  are  insensitive  to  the  differences  between  the  PSF  models  implemented  in  these  IRFs. 
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window  functions  (Eq.  (5))  for  the  PSF  implemented  in  the 
P6_V3  IRFs  used  in  this  analysis  to  those  for  the  PSF  in  the 
more  recently  updated  P6_vil  IRFs.  The  P6_vil  IRFs  use  a 
modified  functional  form  for  the  PSF,  and  for  energies 
above  1  GeV  the  PSF  implemented  in  P6_vil  was  cali¬ 
brated  using  in-flight  data,  while  in  P6_V3  the  PSF  was 
based  on  Monte  Carlo  simulations.  Figure  9  shows  the 
beam  window  functions  for  the  PSF  associated  with  the 
front-  and  back-converting  events  for  each  set  of  IRFs,  at 
the  log-center  of  each  energy  bin  used  in  this  analysis.  The 
small  variation  between  the  window  functions  of  the  two 
IRFs  confirms  that  differences  between  the  PSF  models  in 
these  two  IRFs  are  not  large  enough  to  affect  the  anisotropy 
measurement  on  the  angular  scales  to  which  this  analysis  is 
sensitive. 

E.  Dependence  on  the  latitude  mask 

In  this  analysis  we  apply  a  generous  latitude  mask  to 
reduce  contamination  of  the  data  by  Galactic  diffuse  emis¬ 
sion.  The  mask  is  intended  to  remove  enough  contamina¬ 


tion  so  that  the  measured  angular  power  can  be  attributed  to 
sources  whose  distribution  is  statistically  isotropic  in  the 
sky  region  we  consider,  i.e.,  a  distribution  which  does  not 
show  any  preferred  direction  on  the  sky.  In  particular,  we 
wish  to  exclude  sources  whose  angular  distribution  exhib¬ 
its  a  strong  gradient  with  Galactic  latitude.  The  effective¬ 
ness  of  the  mask  at  reducing  the  contribution  to  the  angular 
power  from  a  strongly  latitude-dependent  component  can 
be  evaluated  by  considering  the  angular  power  spectrum  of 
the  data  as  a  function  of  latitude  cut.  The  results  are  shown 
in  Figs.  10  and  11. 

At  low  multipoles  (€  ;£  100),  increasing  the  latitude  cut 
significantly  reduces  the  angular  power,  indicating  that  in 
this  multipole  range  the  contamination  by  a  strongly 
latitude-dependent  component,  such  as  Galactic  diffuse 
emission,  is  considerable.  For  155  ^  f?  <  254  at  1-2  GeV 
and  2-5  GeV,  the  angular  power  measured  using  the  30° 
latitude  mask  is  noticeably  smaller  than  when  using  the  20° 
latitude  mask.  However,  at  all  energies  there  are  no  signifi¬ 
cant  differences  in  the  angular  power  measured  for  €  >155 
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FIG.  10  (color  online).  Intensity  angular  power  spectra  of  the  data  calculated  with  different  latitude  cuts.  The  point  source  mask  was 
applied  in  addition  to  the  latitude  mask  in  all  cases.  The  differences  between  the  results  masking  \b\  <  30°  (the  default  latitude  cut) 
and  \b\  <  40°  are  small  for  (,  >  155  for  all  four  energy  bins,  demonstrating  that  the  power  observed  in  the  data  at  these  multipoles  is 
not  strongly  correlated  with  a  component  that  has  a  strong  latitude  dependence  in  the  range  30°  <  \b\  <  40° ,  such  as  the  Galactic 
diffuse  emission.  At  energies  above  5  GeV  convergence  is  seen  for  multipoles  €  >  155  even  when  masking  only  \b\  <  20°.  Points 
from  different  data  sets  are  offset  slightly  in  multipole  for  clarity.  The  lowest  multipole  data  point  for  the  \b\  <  20°  mask  in  each  panel 
is  above  the  range  shown  in  the  figure.  Expanded  versions  of  the  top  panels  are  shown  in  Fig.  1 1 . 
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FIG.  1 1  (color  online).  Expanded  versions  of  top  panels  of  Fig.  10. 


using  the  30°  and  40°  latitude  masks,  and  for  energies 
greater  than  5  GeV  the  20°  latitude  mask  also  yields  con¬ 
sistent  results.  We  conclude  that  applying  the  30°  latitude 
mask  is  sufficient  to  ensure  that  no  significant  amount  of  the 
measured  angular  power  at  £  >  155  originates  from  the 
Galactic  diffuse  emission  or  from  any  source  class  that 
varies  greatly  in  the  region  30°  <  \b\  <  40°. 


F.  Effects  of  masking  on  the  power  spectrum 

To  verify  that  the  results  do  not  depend  sensitively  on  the 
angular  radius  of  the  source  mask,  in  Figs.  12  and  13  we 
compare  the  results  when  masking  a  1°  angular  radius 
around  each  source  with  those  when  masking  the  2°  radius 
used  as  the  default  in  this  work. 
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FIG.  12  (color  online).  Intensity  angular  power  spectra  of  the  data  calculated  with  a  mask  excluding  a  1  °  or  2°  angular  radius  around 
each  source;  excluding  a  2°  angular  radius  is  the  default  in  this  analysis.  The  default  latitude  mask  excluding  \b\  <  30°  was  applied  in 
addition  to  the  source  mask  in  all  cases.  At  all  energies  the  angular  power  spectra  obtained  using  the  different  source  mask  radii  are 
consistent  at  €  >  155  (the  multipole  range  of  interest),  and  above  2  GeV  the  results  are  consistent  at  f  >  55.  Expanded  versions  of  the 
top  panels  are  shown  in  Fig.  13. 
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FIG.  13  (color  online).  Expanded  versions  of  top  panels  of  Fig.  12. 


In  the  1-2  GeV  energy  bin  the  results  show  significant 
differences  at  €  <  155,  however  for  €  >  155  (the  multipole 
range  of  interest)  the  angular  power  spectra  for  the  1 0  and 
2°  source  mask  cases  agree  within  the  error  bars.  In  the 
higher  energy  bins  the  angular  power  spectra  in  all  except 
the  first  multipole  bin  (5  ^  f  <  55,  well  below  the  range  of 
interest)  agree  within  the  error  bars.  Since  varying  the 
angular  size  of  the  region  masked  around  each  source 
does  not  significantly  change  the  measured  angular  power 
at  €  >  155,  we  conclude  that  any  features  that  may  be 
induced  in  the  angular  power  spectra  by  the  morphology  of 
the  source  mask  are  confined  to  low  multipoles  and  there¬ 
fore  do  not  affect  the  measurements  of  CP  reported  in  this 
work. 

In  addition,  we  have  confirmed  that  the  angular  power 
spectra  of  the  front-  and  back-converting  events  are  in  good 
agreement  within  each  energy  bin  in  the  multipole  range  of 
interest  (F  >  155),  and  are  generally  consistent  at  €  <  155 
even  in  the  1-2  GeV  energy  bin  where  the  95%  contain¬ 
ment  radius  of  the  PSF  of  the  back-converting  events  is 
comparable  to  the  angular  radius  used  for  the  source  mask. 
Consequently,  although  the  PSF  associated  with  the  back- 
converting  events  is  larger  than  that  of  the  front-converting 
events,  the  consistency  of  their  angular  power  spectra 
implies  that  the  source  masking  is  sufficiently  effective 
even  at  low  energies. 

The  sharp  latitude  cut  used  in  this  analysis  also  has  the 
potential  to  induce  features  in  the  angular  power  spectrum, 
although  these  would  be  expected  to  appear  on  the  large 
angular  scales  characteristic  of  the  morphology  of  the 
mask.  We  therefore  note  that  the  stability  of  the  angular 
power  spectra  at  €  >  155  for  latitude  cuts  masking  at  least 
\b\  <  30°,  discussed  in  Sec.  VIE  and  demonstrated  in 
Figs.  10  and  11,  indicates  that  the  latitude  mask  does  not 
induce  features  in  the  power  spectrum  at  the  angular  scales 
of  interest. 

The  analysis  of  the  simulated  isotropic  component,  pre¬ 
sented  in  Sec.  VI H,  provides  another  means  of  assessing 
the  impact  of  the  mask  on  the  angular  power  spectra.  Since 
the  isotropic  component  should  only  contribute  to  the 


monopole  (€  =  0)  term  of  the  power  spectrum,  statistically 
significant  deviations  from  zero  power  at  €  >  0  can  be 
attributed  to  the  use  of  the  mask.  We  emphasize  that  the 
consistency  of  the  angular  power  of  the  isotropic  compo¬ 
nent  with  zero  at  €  >  155  indicates  that,  despite  the  com¬ 
plex  morphology  of  the  total  mask,  the  mask  does  not 
induce  features  in  the  angular  power  spectrum  at  the  multi¬ 
poles  of  interest  (f  >  155). 

G.  Dependence  on  the  set  of  masked  sources 

The  recently-released  second  Fermi  LAT  source  catalog 
(2FGL)  [62]  is  an  update  to  the  1FGL  catalog  used  to 
define  the  default  source  mask  adopted  in  this  work.  The 
2FGL  catalog  reports  the  detection  of  1873  sources,  com¬ 
pared  to  the  1451  included  in  the  1FGL  catalog. 

We  briefly  comment  that  one  motivation  for  using  the 
1FGL  catalog,  rather  than  the  2FGL  catalog,  to  define  the 
source  mask  in  our  default  analysis  is  that  the  1FGL 
catalog  was  also  used  in  the  Fermi  LAT  source  count 
distribution  analysis  [10].  The  results  of  that  study  are 
closely  related  to  the  interpretation  of  the  results  of  the 
current  analysis,  and  so  our  choice  to  mask  that  same 
source  list  in  our  default  analysis  allows  the  results  of  the 
two  analyses  to  be  used  together  straightforwardly. 
However,  it  is  natural  to  ask  to  what  extent  the  measured 
angular  power  reported  in  the  data  may  be  attributable  to 
the  additional  sources  resolved  in  the  2FGL  catalog. 

We  address  this  question  by  analyzing  the  data  using  a 
source  mask  defined  by  the  2FGL  sources  and  comparing 
the  results  to  those  obtained  using  the  1FGL  source  mask. 
We  repeat  the  analysis  of  the  data  using  the  default  pipe¬ 
line,  changing  only  the  source  mask;  the  total  mask  is 
defined  by  the  source  mask  combined  with  the  default 
latitude  cut  masking  \b\  <  30°.  When  combined  with  the 
default  latitude  cut,  the  2FGL  source  mask  results  in  an 
unmasked  sky  fraction  /sky  =  0.295,  a  small  decrease 
compared  to  /sky  =  0.325  when  using  the  1FGL  source 
mask. 

The  angular  power  spectra  of  the  data  analyzed  using  the 
2FGL  catalog  to  define  the  source  mask  are  shown  in 
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FIG.  14  (color  online).  Intensity  angular  power  spectra  of  the  data  calculated  using  the  source  mask  defined  by  the  2FGL  catalog 
compared  with  the  results  using  the  1FGL  catalog;  the  source  mask  defined  by  the  1FGL  catalog  is  the  default  used  in  this  analysis. 
The  angular  power  at  €  >  155  is  smaller  in  the  2FGL  case  by  —20-30%  in  the  bins  spanning  1-10  GeV  and  by  ~70%  at  10-50  GeV. 
Expanded  versions  of  the  top  panels  are  shown  in  Fig.  15. 


Figs.  14  and  15,  compared  with  the  results  of  the  default 
data  analysis  which  uses  the  1FGL  catalog.  The  angular 
power  CP  measured  in  the  data  using  the  2FGL  source 
mask  is  reduced  relative  to  the  1FGL  case  (see  Table  I), 
while  the  measurement  uncertainties  remain  roughly  the 
same  as  in  the  1FGL  case.  The  decrease  in  CP  is  —20-30% 
in  the  1-2,  2-5,  and  5-10  GeV  energy  bins,  however 


significant  detections  (  >  3cr)  are  still  found  in  these  three 
bins.  A  —  70%  decrease  in  CP  is  seen  in  the  10-50  GeV 
bin,  and  due  to  the  large  measurement  uncertainty  the 
significance  of  the  measurement  in  this  bin  falls  from 
2.7o-  to  0.8(7.  The  significance  of  the  detected  fluctuation 
angular  power  over  all  four  energy  bins  remains  greater 
than  7cr. 
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FIG.  15  (color  online).  Expanded  versions  of  top  panels  of  Fig.  14. 
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We  can  estimate  the  expected  decrease  in  angular 
power  when  masking  the  2FGL  sources  by  calculating 
the  difference  in  angular  power  produced  when  the  source 
detection  threshold  is  reduced  from  the  1FGL  to  the  2FGL 
catalog  level,  following  the  approach  used  to  calculate  the 
angular  power  of  the  simulated  sources  in  Sec.  VI B.  We 
assume  the  sources  follow  the  flux  distribution  function 
dN/dS  given  by  Eq.  (10)  with  the  same  parameters  given 
in  that  section  as  the  best-fit  for  the  high-latitude 
Fermi  sources  in  the  1.04-10.4  GeV  band.  Calculating 
CP  via  Eq.  (11)  for  an  assumed  flux  threshold  of  ~5  X 
10-10  cm“2s-1,  appropriate  for  the  1FGL  catalog  [50], 
yields  CP~9.4X  lCT18(cm“2s_1sr_1GeV_1)2sr.  Using  a 
lower  flux  threshold  of  ~4  X  10“ 10  cm-2  s-1,  appropriate 
for  the  2FGL  catalog  [62],  gives  CP~6.8X10-18 
(cm_2s_1sr^1GeV_1)2sr7  which  is  indeed  a  roughly  30% 
decrease  in  CP,  as  observed  in  the  data. 

H.  Comparison  of  data  and  simulated  models 

To  understand  the  origin  of  the  angular  power  measured 
in  the  data,  we  compare  the  angular  power  spectra  of  the 
default  data  to  those  of  the  default  simulated  model  and  the 
alternate  simulated  model,  described  in  Sec.  V.  The  simu¬ 


lated  models  were  processed  and  their  angular  power 
spectra  calculated  using  the  same  analysis  pipeline  as  the 
data,  and  thus  we  expect  the  angular  power  spectra  of  the 
data  and  models  to  be  consistent  if  the  models  accurately 
reflected  the  statistical  properties  of  the  emission  on  the 
relevant  angular  scales. 

Figures  16  and  17  present  the  angular  power  spectra  of 
the  data  and  models.  The  angular  power  spectra  of  the 
two  models  agree  very  well  at  all  energies  at  multipoles 
above  €  =  105.  At  all  energies  and  scales,  both  models 
exhibit  less  angular  power  than  the  data.  Moreover,  the 
amplitude  of  the  detected  angular  power  in  both  models  is 
inconsistent  with  that  of  the  data  at  >95%  CL  in  the  three 
energy  bins  spanning  1-10  GeV,  and  at  >90%  CL  in  the 
10-50  GeV  bin  (see  Table  III).  The  lack  of  significant 
power  at  high  multipoles  in  either  simulated  model  indi¬ 
cates  that  the  Galactic  diffuse  emission,  as  implemented  in 
these  models,  provides  a  negligible  contribution  to  the 
anisotropy  at  €  >  155.  At  lower  multipoles,  the  discrep¬ 
ancy  between  the  data  and  models  and  between  the  two 
models  may  be  due  to  the  presence  of  large-scale  features 
in  the  data  which  are  not  included  in  the  models,  however 
we  defer  a  full  investigation  of  the  origin  of  the  low- 
multipole  angular  power  to  future  work. 
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FIG.  16  (color  online).  Angular  power  spectra  of  the  data,  the  default  simulated  model  (MODEL),  and  the  alternate  simulated  model 
(ALT  MODEL).  The  angular  power  spectra  of  the  two  models  are  in  good  agreement  in  all  energy  bins.  The  smaller  amplitude  angular 
power  at  €  >  155  measured  at  lower  significance  in  both  models  is  inconsistent  with  the  angular  power  observed  in  the  data  at  all  energies. 
Points  from  different  data  sets  are  offset  slightly  in  multipole  for  clarity.  Expanded  versions  of  the  top  panels  are  shown  in  Fig.  17. 
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FIG.  17  (color  online).  Expanded  versions  of  top  panels  of  Fig.  16. 


The  contributions  to  the  angular  power  spectrum  of  the 
individual  components  of  the  default  model  are  shown  in 
Figs.  18  and  19.  At  all  energies  the  only  component  con¬ 
tributing  significantly  to  the  total  power  is  the  Galactic 
diffuse  emission.  The  contribution  from  the  isotropic  com¬ 
ponent  is  negligible,  since  this  component  is  isotropic  by 
construction  and  thus,  after  the  photon  noise  is  subtracted, 
it  should  only  contribute  to  the  monopole  (f  =  0)  term. 
The  deviations  from  zero  of  the  isotropic  component  in  the 
lowest  multipole  bin  (5  ^  F  <  54)  may  be  due  to  imper¬ 
fect  correction  of  the  effects  of  the  mask  in  this  multipole 
regime.  The  source  catalog  component  contributes  zero 
power  at  all  energies  and  multipoles  since  the  emission 
maps  of  this  simulated  component  contain  only  events 
from  sources  which  are  masked  in  the  analysis.  The  con¬ 
sistency  of  the  source  catalog  angular  power  with  zero 
indicates  that  the  source  masking  is  effective.  We  remark 
that  in  general  the  angular  power  spectra  of  distinct  com¬ 
ponents  are  not  linearly  additive  due  to  contributions  from 
cross-correlations  between  the  components.  The  total 
power  of  the  model  is,  however,  very  consistent  with  the 
total  power  in  the  Galactic  component,  with  slight  discrep¬ 
ancies  likely  arising  from  masking  effects,  since  the 
Galactic  and  isotropic  components  should  have  no  cross¬ 
correlation  power  and  the  simulated  sources  were  fully 
masked. 


TABLE  III.  Significance  of  the  difference  ACP  between  inten¬ 
sity  angular  power  CP  for  155  ^  ^  504  in  the  default  data  and 

the  default  simulated  model  in  each  energy  bin.  The  associated 
measurement  uncertainties  can  be  taken  to  be  Gaussian. 


F 

^min 

F 

^max 

Significance  of  ACP 

1.04 

1.99 

3.5<r 

1.99 

5.00 

4.5cr 

5.00 

10.4 

2.  Ocr 

10.4 

50.0 

1.7cr 

VII.  ENERGY  DEPENDENCE  OF  ANISOTROPY 
IN  THE  DATA 

The  energy  dependence  of  the  fluctuation  angular  power 
can  be  used  to  identify  the  presence  of  multiple  distinct 
contributors  to  the  emission  [49].  Because  the  fluctuation 
angular  power  characterizes  only  the  angular  distribution 
of  the  emission,  independent  of  the  intensity  normaliza¬ 
tion,  it  is  exactly  energy-independent  for  a  single  source 
class  as  long  as  the  members  of  the  class  have  the  same 
observed  energy  spectrum.  In  general,  the  fluctuation 
angular  power  of  a  single  source  class  may  show  energy 
dependence  due  to  large  variation  of  the  energy  spectra  of 
individual  sources  within  a  population,  and,  for  cosmologi¬ 
cal  source  classes,  the  effects  of  redshifting  and  attenuation 
of  high-energy  gamma  rays  by  the  extragalactic  back¬ 
ground  light  (EBL).  Redshifting  and  EBL  attenuation  are 
expected  to  be  important  only  for  populations  for  which  a 
significant  fraction  of  the  observed  intensity  originates 
from  high-redshift  members,  with  EBL  attenuation  rele¬ 
vant  only  at  observed  energies  of  several  tens  of  GeV.  All 
of  these  effects  are  most  prominent  when  the  source  spec¬ 
tra  have  hard  features  such  as  lines  or  cut-offs;  smoothly- 
varying  source  spectra,  such  as  power-law  energy  spectra, 
typically  generate  more  mild  energy  dependence  in  the 
fluctuation  angular  power. 

The  fluctuation  anisotropy  energy  spectrum  of  the  data 
is  shown  in  the  top  panel  of  Fig.  20.  The  fluctuation 
angular  power  CP/(/)2  in  each  energy  bin  was  obtained 
by  weighted  averaging  of  the  unbinned  fluctuation 
angular  power  spectrum  over  155  <  F  <  504,  weighting 
the  measured  angular  power  at  each  multipole  by  its  mea¬ 
surement  uncertainty;  these  values  are  reported  in  Table  I. 
Each  point  is  located  at  the  logarithmic  center  of  the  energy 
bin. 

A  power-law  fit  of  the  fluctuation  angular  power  as  a 
function  of  energy  CP/(/)2  «  E~Tv  yields  TF  =  0.076  ± 
0.139  (—0.082  ±  0.158  for  the  cleaned  data),  consistent 
with  no  energy  dependence  over  the  energy  range  consid- 
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FIG.  18  (color  online).  Angular  power  spectra  ol  the  components  of  the  default  simulated  model  (MODEL).  As  expected,  most  of 
the  total  angular  power  at  all  multipoles  (TOTAL  MODEL)  is  due  to  the  GAL  component.  By  construction,  the  isotropic  component 
(ISO)  component  contributes  no  significant  angular  power;  likewise,  the  source  component  (CAT)  provides  no  contribution  because  all 
sources  were  masked.  Points  corresponding  to  the  TOTAL  MODEL  are  offset  slightly  in  multipole  for  clarity.  Expanded  versions  of 


the  top  panels  are  shown  in  Fig.  19. 

ered.  The  best-fit  constant  value  of  CP/(/)2  across  all  four 
energy  bins  is  9.05 ± 0.84 X  10_6sr  (6.94±0.84X  10_6sr 
for  the  cleaned  data).  The  results  of  these  fits  for  the  data 
with  and  without  foreground  cleaning  are  summarized  in 
Table  IV,  along  with  the  results  for  the  energy  dependence 
of  the  intensity  angular  power,  discussed  below.  The  lack 
of  a  clear  energy  dependence  in  the  fluctuation  angular 
power  is  consistent  with  a  single  source  class  providing  the 


dominant  contribution  to  the  anisotropy  and  a  constant 
fractional  contribution  to  the  intensity  over  the  energy 
range  considered,  although  due  to  the  large  measurement 
uncertainties  contributions  from  additional  source  classes 
cannot  be  excluded.  This  is  especially  true  for  sources 
whose  contribution  to  the  intensity  peaks  at  £  S 
10  GeV.  Furthermore,  due  to  the  coarseness  of  the  energy 
binning,  this  analysis  is  not  sensitive  to  features  in  the 
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FIG.  19  (color  online).  Expanded  versions  of  top  panels  of  Fig.  18. 


083007-22 


ANISOTROPIES  IN  THE  DIFFUSE  GAMMA-RAY  . . . 
2.0 1CT5 

1.5*1 0-5 

F  1.O10-5 

A  5.0-Kr6 

0 

-5.0KT6 

1  10 
Energy  [GeV] 


DATA  X 
DATA:CLEANED  □ 


> 

tu 

O 


10" 


10"s  - 


10"“  - 


S'  to-22  - 
<1 


10- 


DATA  X 
DATA:CLEANED  □ 


10 

Energy  [GeV] 


FIG.  20  (color  online).  Anisotropy  energy  spectra  of  the  data. 
Top:  Fluctuation  anisotropy  energy  spectrum.  The  data  are 
consistent  with  no  energy  dependence  over  the  energy  range 
considered,  although  a  mild  energy  dependence  is  not  excluded. 
Bottom:  Differential  intensity  anisotropy  energy  spectrum.  The 
energy  dependence  is  consistent  with  that  arising  from  a  single 
source  population  with  a  power-law  intensity  energy  spectrum 
with  spectral  index  Ts  =  2.40  ±  0.07  for  the  default  data 
(2.33  ±  0.08  for  the  cleaned  data). 


anisotropy  energy  spectrum  localized  to  narrow  energy 
bands. 

If  a  single  source  class  dominates  the  anisotropy  at  all 
energies  considered,  the  differential  intensity  angular 
power  spectrum  Cf/(AE)2  scales  with  energy  as  the  inten¬ 
sity  energy  spectrum  squared  ( dN/dE )2  of  that  source 
class.  For  example,  for  a  source  class  with  a  power-law 
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photon  spectrum  dN/dE  E~r\  Q/(A£)2  a  E~2T\ 
We  can  therefore  use  this  energy  scaling  to  constrain 
the  energy  spectrum  of  the  dominant  contributor  to  the 
anisotropy,  under  the  assumption  that  the  measured  angular 
power  (but  not  necessarily  the  total  measured  intensity) 
originates  from  a  single  source  class. 

Here  we  obtain  the  differential  intensity  angular  power 
CP/(A£’)2  by  dividing  the  intensity  angular  power  CP  in 
each  energy  bin  by  the  bin  size  squared.  The  differential 
intensity  anisotropy  energy  spectrum  of  the  data  is  shown 
in  the  bottom  panel  of  Fig.  20.  The  CP  are  the  best-fit 
values  for  155  <  L  <  504,  i.e.,  the  weighted  average  of 
Cf  in  that  multipole  range,  reported  in  Table  I,  and  each 
data  point  is  located  at  the  logarithmic  center  of  the  energy 
bin.  The  results  of  fitting  CP/(A£)2  E  1 1  are  given  in 
Table  IV.  Identifying  I  ,  =  2TS,  the  best  fit  of  the  energy 
dependence  suggests  that  the  anisotropy  is  contributed  by  a 
source  class  with  a  power-law  photon  spectrum  character¬ 
ized  by  Ts  =  2.40  ±  0.07  (2.33  ±  0.08  for  the  cleaned 
data),  assuming  only  one  source  class  contributes  appreci¬ 
ably  to  the  anisotropy.  As  the  single  power-law  energy 
dependence  provides  a  very  good  fit  to  the  data,  attributing 
the  anisotropy  to  a  single  source  class  is  a  plausible 
interpretation. 

We  note  that  the  spectral  index  implied  for  the  dominant 
source  class  contributing  to  the  anisotropy  is  in  excellent 
agreement  with  the  mean  intrinsic  spectral  index  of  blazars 
as  inferred  from  the  Fermi-detected  members  [10], 
strongly  supporting  the  interpretation  of  the  measured 
anisotropy  as  originating  from  unresolved  blazars.  We 
caution,  however,  that  due  to  the  variation  between  indi¬ 
vidual  blazars’  spectral  indices,  as  well  as  possible  effects 
of  EBL  attenuation  and  redshifting,  the  fluctuation  angular 
power  from  blazars  could  exhibit  some  energy  dependence 
in  the  range  considered  here.  Therefore,  assuming  that 
blazars  are  the  dominant  source  class  contributing  the 
anisotropy  could  lead  to  tension  with  the  flatness  of  the 
measured  fluctuation  anisotropy  energy  spectrum. 
Additional  support  for  a  blazar  interpretation  could  be 
provided  by  a  detailed  study  of  the  energy-dependent 
anisotropy  arising  from  specific  blazar  population  models, 
calibrated  to  match  the  properties  of  Fermi-detected  blaz¬ 
ars,  and  the  consistency  of  the  predicted  anisotropy  of  these 
models  with  the  measured  amplitude  of  the  angular  power. 
We  defer  a  careful  treatment  of  this  subject  to  future  work. 


TABLE  IV.  Energy  dependence  of  angular  power  for  155  <  €  <  504  in  each  energy  bin  for  the  data  processed  with  the  default 
analysis  pipeline  and  the  Galactic-foreground-cleaned  data.  The  best-fit  constant  value  of  the  fluctuation  angular  power  (CP/(/)2)  over 
1-50  GeV  is  obtained  by  weighted  averaging  of  CP/(/)2  of  the  four  energy  bins.  The  best-fit  parameters  and  associated  x~  per  degree 
of  freedom  (d.o.f.)  are  given  for  fits  of  the  fluctuation  angular  power  to  CP/(/}2  =  AF(E / E0)~rf  and  the  differential  intensity  angular 
power  to  CP/(A£)2  =  Ai(E / E0)~Ti ,  with  E0  =  1  GeV.  The  value  of  Ai  is  given  in  terms  of  AT/AT0,  where  AI0  =  10-18 
(cm 2s 1  sr"1  GeV-1)2  sr. 


<CP/</>2)  [10"6  sr] 

AF  [10"6  sr] 

rF 

*2/d.o.f. 

Ai/AI0 

r, 

;r2/d.o.f. 

DATA 

9.05  ±  0.84 

9.85  ±  1.73 

0.076  ±  0.139 

0.41 

45.1  ±  7.8 

4.79  ±  0.13 

0.19 

DATA :  CLEANED 

6.94  ±  0.84 

6.31  ±  1.44 

-0.082  ±  0.158 

0.12 

29.4  ±  6.6 

4.66  ±0.15 

0.035 
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VIII.  DISCUSSION 

Prior  work  has  generated  predictions  for  the  angular 
power  spectra  of  several  source  populations  which  may 
contribute  to  the  IGRB.  In  most  cases,  the  predictions  for 
the  anisotropy  of  the  emission  from  a  single  source  class 
have  been  cast  in  terms  of  fluctuation  angular  power 
Cf  /(I)2,  where  Cf  is  the  intensity  angular  power  spectrum 
of  the  source  class  and  (I)  its  mean  collective  intensity  in  a 
specified  energy  range.  Since  the  intensity  contributions  of 
most  gamma-ray  source  classes  to  the  IGRB  are  subject  to 
large  uncertainties,  it  is  convenient  to  consider  the  fluctua¬ 
tion  angular  power,  since  this  quantity  is  independent  of 
the  overall  normalization  of  the  intensity.  This  convention 
is  particularly  useful  when  the  spatial  distribution,  number 
density,  and  relative  flux  distribution  of  the  sources  is 
known  or  modeled,  and  the  uncertainty  in  the  collective 
intensity  can  be  translated  into  a  multiplicative  factor  that 
uniformly  scales  the  observed  intensity  in  all  sky  direc¬ 
tions.  For  this  reason,  the  fluctuation  angular  power  is  very 
well  suited  for  characterizing  an  indirect  dark  matter  signal 
since  the  intensity  normalization  scales  linearly  with  the 
assumed  annihilation  cross-section  or  decay  rate. 

By  comparing  the  measured  fluctuation  angular  power 
with  predictions  for  various  source  classes,  we  can  place 
constraints  on  the  fractional  contribution  from  each  source 
class  to  the  total  intensity  by  requiring  that  the  fluctuation 
angular  power  of  the  total  emission  is  not  exceeded. 
Assuming  that  each  contributing  source  class  is  uncorre¬ 
lated  with  the  others  and  the  Poisson  component  dominates 
the  angular  power  spectrum  of  each  source  class  at  the 
multipoles  considered,  the  intensity  angular  power  of  the 
total  emission  is  given  by 

Cp.tot  =  Up  i  +  Cp2  +  . . .  (12) 


and  so  the  fluctuation  angular  power  of  the  total  intensity  is 


Up, tot  _  Up  1  Up  ^ 

(hot)2  (hot)2  (hot)2 


(13) 


Rewriting  the  fractional  contribution  from  source  class  i  to 
the  total  intensity  /,  =  (h)/(hot), 


Up,  tot 

<4^ 


r2  UP,1  , 

~  w 
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2/r  \2  + 
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(14) 


If  we  allow  a  single  source  class  i  to  contribute  all  of  the 
measured  angular  power,  the  source  class  is  constrained 
such  that 


Up, tot/ (hot)2 

cP,,/(h)2 


(15) 


Source  classes  whose  predicted  fluctuation  angular  power 
exceeds  the  measured  fluctuation  angular  power  therefore 
cannot  contribute  the  entirety  of  the  measured  intensity. 

It  is  important  to  note  that  the  total  intensity  of  the  IGRB 
in  this  analysis  is  not  equivalent  to  the  intensity  of  the 
isotropic  emission  reported  in  [5],  since  that  analysis 
employed  much  more  stringent  selection  cuts  to  remove 
charged  particle  contamination,  and  used  a  fitting  proce¬ 
dure  to  remove  contributions  from  nonisotropic  compo¬ 
nents.  In  this  analysis  the  total  intensity  of  the  emission  is 
simply  the  intensity  that  remains  after  the  mask  is  applied, 
which  may  include  some  emission  from  nonisotropic  com¬ 
ponents,  as  well  as  a  non-negligible  amount  of  charged 
particle  contamination.  However,  we  emphasize  that  since 
the  charged  particle  contamination  is  presumed  to  be 
nearly  isotropic,  with  any  potential  fluctuations  confined 
to  large  angular  scales,  it  should  not  contribute  to  the 
intensity  angular  power  at  the  multipoles  considered  here 
(f  >  155),  and  so  a  more  robust  comparison  of  models 
with  the  data  could  be  achieved  by  comparing  the  predicted 
intensity  angular  power  to  the  measurement. 

We  now  compare  our  measurement  to  existing  predic¬ 
tions  from  the  literature  for  the  angular  power  spectra  of 
various  gamma-ray  source  classes,  and  summarize  these 
results  in  Table  V.  We  caution  that  the  predicted  angular 
power  can  depend  sensitively  upon  the  adopted  source 
model  (in  particular  the  shape  of  the  flux  distribution), 
the  assumed  source  detection  threshold,  and,  for  cosmo¬ 
logical  source  classes,  assumptions  regarding  the  effect 
on  the  observed  energy  spectrum  of  attenuation  of 
high-energy  photons  by  interactions  with  the  EBL. 
Consequently,  the  constraints  derived  in  this  section  should 


TABLE  V.  Maximum  fractional  contribution  of  various  source  populations  to  the  IGRB  intensity  that  is  compatible  with  the  best-fit 
constant  value  of  the  measured  fluctuation  angular  power  in  all  energy  bins,  (CP/(I)2)  =  9.05  X  10-6  sr  for  the  default  data  analysis  or 
(CP/(/)2}  =  6.94  X  10-6  sr  for  the  Galactic-foreground-cleaned  data  analysis.  Indicative  values  for  the  fluctuation  angular  power 
Cf/(I)2  of  each  source  class  are  taken  from  existing  literature  (see  text  for  details)  and  evaluated  at  f  =  100. 


Source  class 

Predicted  C\qq/{I)2 

Maximum  fraction  of  IGRB  intensity 

[sr] 

DATA 

DATA :  CLEANED 

Blazars 

2  X  10“4 

21% 

19% 

Star-forming  galaxies 

2  X  10“7 

100% 

100% 

Extragalactic  dark  matter  annihilation 

1  X  10“5 

95% 

83% 

Galactic  dark  matter  annihilation 

5  X  10“5 

43% 

37% 

Millisecond  pulsars 

3  X  10“2 

1.7% 

1.5% 
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be  taken  only  as  indicative  values  for  these  source 
populations. 

Reference  [27]  predicted  the  fluctuation  anisotropy  from 
unresolved  blazars  CP/(/)2  ~  2  X  10-4  sr  at  €  ~  100  (see 
Fig.  4  of  that  work).  This  is  a  factor  of  —20  larger  than  the 
fluctuation  angular  power  of  — 10-5  sr  measured  in  the 
data,  which  suggests  that  emission  from  blazars,  assuming 
the  model  adopted  in  that  study,  contributes  less  than 
—  1/ V20  —  20%  of  the  total  intensity.  Note,  however, 
that  the  flux  threshold  for  sources  in  the  1FGL  catalog  is 
between  0.5  and  1  X  10-9  photons  cm“2  s-1  for  \b\  >  30°, 
higher  than  the  threshold  assumed  in  [27].  If  the  blazar 
luminosity  function  is  identical  to  the  one  assumed  in  [27], 
this  discrepancy  in  thresholds  would  imply  that  the  pre¬ 
diction  for  the  blazar  anisotropy  in  [27]  is  underestimated 
with  respect  to  the  one  applicable  to  our  analysis,  since  our 
masked  maps  include  more  bright  unresolved  blazars.  As  a 
result,  the  constraint  on  the  fractional  intensity  contribu¬ 
tion  to  the  IGRB  from  blazars  for  this  model  from  our 
measurement  would,  if  anything,  be  stronger. 

In  contrast  to  the  larger  anisotropy  expected  from  blaz¬ 
ars,  the  fluctuation  angular  power  at  £  —  100  predicted  for 
star-forming  galaxies  by  Ref.  [30]  is  —2  X  10~7  sr  at 
1  GeV,  far  below  the  value  measured  in  this  analysis. 
Since  star-forming  galaxies  would  thus  provide  a  subdo¬ 
minant  contribution  to  the  measured  angular  power,  this 
anisotropy  measurement  does  not  constrain  their  contribu¬ 
tion  to  the  total  IGRB  intensity. 

The  anisotropy  from  dark  matter  annihilation  in  extra- 
galactic  structures  is  predicted  to  be  slightly  smaller  than 
that  from  unresolved  blazars,  although  estimates  can  vary 
substantially  due  to  differences  in  the  adopted  models. 
Moreover,  for  extragalactic  dark  matter  annihilation  the 
amplitude  of  the  expected  anisotropy  can  be  highly  sensi¬ 
tive  to  the  energy  spectrum  of  the  emission.  The  source 
energy  spectrum  depends  on  the  dark  matter  particle  mass 
and  dominant  annihilation  channels,  while  the  observed 
energy  spectrum  is  affected  by  redshifting  and  EBL  at¬ 
tenuation.  These  factors  can  introduce  a  nontrivial  energy 
dependence  into  the  amplitude  of  the  anisotropy,  particu¬ 
larly  for  high  mass  (  —  1  TeV)  dark  matter  candidates.  As 
a  benchmark  range,  Refs.  [23,27,39]  predict  the  anisotropy 
from  annihilation  of  extragalactic  dark  matter  to  be 
— 10— 6— 10— 5  sr  at  €  —  100  at  energies  of  a  few  GeV, 
comparable  to  the  measured  value. 

The  anisotropy  from  annihilation  in  Galactic  dark  matter 
substructure  is  expected  to  be  much  larger  than  that  from 
extragalactic  dark  matter.  While  variations  in  the  assumed 
properties  of  Galactic  substructure  can  lead  to  order-of- 
magnitude  or  larger  variations  in  the  predicted  angular 
power,  for  typical  assumptions  the  predicted  fluctuation 
angular  power  is  —5  X  10-5  sr  at  €  —  100  (e.g..  Model  A1 
in  Ref.  [33]),  which  implies  that  dark  matter  annihilation 
can  contribute  less  than  —43%  of  the  total  intensity. 
However,  adopting  alternative  models  for  the  substructure 
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properties  can  increase  or  decrease  the  predicted  angular 
power  by  as  much  as  —2  orders  of  magnitude  [32-34],  so 
the  measured  angular  power  represents  a  strong  constraint 
on  some  substructure  models. 

Galactic  gamma-ray  millisecond  pulsars  (MSPs)  have 
also  been  considered  as  possible  contributors  to  the  inten¬ 
sity  and  anisotropy  of  the  IGRB  due  to  their  extended 
latitude  distribution  [15,31].  The  emission  from  Galactic 
MSPs  is  expected  to  feature  very  large  fluctuation  anisot¬ 
ropy  due  to  the  relatively  low  number  density  of  this  source 
class  compared  to  dark  matter  substructure  or  extragalactic 
source  populations.  Ref.  [31]  predicts  fluctuation  angular 
power  at  high  Galactic  latitudes  of  —0.03  sr  at  €  —  100  for 
this  Galactic  source  class,  which  implies  a  contribution  to 
the  total  IGRB  intensity  of  no  more  than  a  few  percent. 

In  addition  to  the  specific  source  populations  considered 
in  this  section,  other  Galactic  source  populations  for  which 
anisotropy  predictions  do  not  yet  exist  in  the  literature  may 
also  contribute  to  the  anisotropy  as  well  as  the  intensity  of 
the  high-latitude  diffuse  emission.  These  include  normal 
pulsars,  as  well  as  populations  currently  too  faint  to  have 
had  individual  members  detected  by  Fermi.  The  properties 
of  these  populations  can  be  constrained  by  both  low- 
latitude  and  high-latitude  source  count  analysis  (in  the 
case  that  individual  members  have  been  detected)  [63], 
and  also  by  the  anisotropy  analysis  described  in  this  study. 
We  leave  the  detailed  study  of  this  to  future  work. 

We  note  that  constraints  derived  in  this  section  have  not 
taken  into  account  information  about  the  likely  energy 
spectrum  of  the  dominant  contributing  population,  dis¬ 
cussed  in  Sec.  VII,  which  is  incompatible  with  sources 
known  or  expected  to  feature  spectral  peaks  at  the  energies 
we  consider  (for  example,  Galactic  and  extragalactic  dark 
matter  and  MSPs).  A  careful  study  combining  all  observ¬ 
ables  obtained  in  this  work  would  almost  certainly 
yield  stronger  constraints  on  contributing  populations. 
Furthermore,  we  have  discussed  the  constraints  obtainable 
on  specific  source  populations  by  requiring  that  the  total 
anisotropy  from  each  population  does  not  exceed  the  mea¬ 
sured  value.  We  emphasize,  however,  that  stronger  bounds 
could  be  derived  if  some  fraction  of  the  total  anisotropy 
could  be  robustly  attributed  to  one  or  more  confirmed 
source  classes,  thereby  reducing  the  anisotropy  available 
to  additional  contributors. 

IX.  CONCLUSIONS 

The  statistical  properties  of  the  IGRB  encode  detailed 
information  about  the  origin  of  this  emission.  The  ad¬ 
vanced  capabilities  of  the  Fermi  LAT,  most  notably  its 
improved  angular  resolution  and  large  effective  area, 
have  enabled  a  sensitive  measurement  of  small  angular- 
scale  anisotropies  in  the  IGRB.  Using  —22  months  of  data, 
we  performed  an  angular  power  spectrum  analysis  of  the 
high-latitude  diffuse  emission  measured  by  the  Fermi  LAT. 
Significant  angular  power  above  the  photon  noise  level  is 
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detected  in  the  data  at  multipoles  155  <  F  <  504  in  three 
energy  bins  spanning  1-10  GeV,  and  is  measured  at  lower 
significance  in  the  10-50  GeV  energy  bin.  The  primary 
limitation  of  the  measurement  at  high  energies  is  low  event 
statistics,  which  results  in  the  measurement  uncertainties 
being  dominated  by  the  photon  noise.  In  this  regime  the 
measurement  uncertainties  scale  roughly  inversely  to  the 
number  of  events,  and  hence  increasing  the  statistics  by  a 
factor  of  2  or  3  could  lead  to  a  large  enough  improvement 
in  the  sensitivity  of  the  analysis  to  allow  a  confident 
detection  of  angular  power  in  this  energy  range  and  greater 
sensitivity  to  energy-dependent  anisotropy. 

The  angular  power  measured  in  the  data  at  155  ^  < 

504  is  consistent  with  a  constant  value  within  each  energy 
bin,  and  the  scale  independence  of  the  signal  suggests  that 
it  originates  from  one  or  more  unclustered  populations  of 
point  sources.  Comparing  the  measured  angular  power 
with  predictions  for  known  and  proposed  gamma-ray 
source  classes,  constraints  can  be  obtained  on  the  collec¬ 
tive  intensity  and  properties  of  source  populations  that 
contribute  to  the  IGRB.  The  fluctuation  angular  power 
detected  in  this  analysis  falls  below  the  level  predicted  for 
many  source  classes,  including  blazars,  MSPs,  and  some 
scenarios  for  dark  matter  annihilation  in  Galactic  and  ex- 
tragalactic  structures.  In  these  cases  the  measured  ampli¬ 
tude  of  the  fluctuation  angular  power  limits  the  contribution 
to  the  total  IGRB  intensity  from  each  source  class. 

The  measured  fluctuation  angular  power  is  consistent 
with  a  constant  value  over  the  energy  range  considered, 
however,  due  to  the  relatively  large  measurement  uncer¬ 
tainties  and  limited  number  of  energy  bins,  a  mild  energy 
dependence  in  this  quantity  cannot  be  excluded.  The 
absence  of  a  strong  energy  dependence  in  the  fluctuation 
anisotropy  energy  spectrum  suggests  that  a  single  source 
class  may  provide  the  dominant  contribution  to  the  anisot¬ 
ropy  while  providing  a  constant  fractional  contribution  to 
the  intensity  of  the  IGRB  over  the  energy  range  considered. 
We  caution,  however,  that  this  analysis  is  not  sensitive  to 
structure  in  the  anisotropy  energy  spectrum  that  is  confined 
to  small  energy  ranges,  since  the  requirement  of  large  event 
statistics  to  detect  anisotropies  at  the  measured  level  pre¬ 
cludes  fine  energy  binning  of  the  data.  We  anticipate  that 
future  analyses  that  draw  on  larger  data  sets  will  be  more 
sensitive  to  localized  features  in  the  anisotropy  energy 
spectrum. 

The  energy  dependence  of  the  intensity  angular  power  of 
the  data  is  well-described  by  that  arising  from  a  single 
source  class  with  a  power-law  photon  spectrum  with  index 
Ts  =  2.40  ±  0.07.  Interestingly,  this  value  closely  matches 
the  mean  intrinsic  spectral  index  for  blazars  as  determined 
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from  recent  Fermi  LAT  measurements.  While  alternative 
scenarios  invoking  contributions  from  more  than  one 
source  class  to  explain  the  energy  dependence  of  the 
angular  power  are  in  principle  possible,  the  interpretation 
of  the  measured  power  as  originating  from  a  single  source 
class  with  a  power-law  energy  spectrum  is  an  excellent  fit 
to  the  data.  To  identify  a  specific  population  or  populations 
as  the  source  of  the  measured  IGRB  anisotropy,  detailed 
analysis  of  population  models  for  plausible  source  classes 
will  be  essential  in  order  to  verify  that  both  the  predicted 
intensity  energy  spectrum  of  the  IGRB  and  the  correspond¬ 
ing  anisotropy  signal  provide  a  consistent  explanation  of 
the  data. 
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